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We explore the thermal properties of hot and dense matter using a model that reproduces the 
empirical properties of isospin symmetric and asymmetric bulk nuclear matter, optical model fits 
to nucleon-nucleus scattering data, heavy-ion flow data in the energy range 0.5-2 GeV/A, and 
the largest well-measured neutron star mass of 2 M 0 . Results of this model which incorporates 
finite range interactions through Yukawa type forces are contrasted with those of a zero-range 
Skyrme model that yields nearly identical zero-temperature properties at all densities for symmetric 
and asymmetric nucleonic matter and the maximum neutron star mass, but fails to account for 
heavy-ion flow data due to the lack of an appropriate momentum dependence in its mean field. 
Similarities and differences in the thermal state variables and the specific heats between the two 
models are highlighted. Checks of our exact numerical calculations are performed from formulas 
derived in the strongly degenerate and non-degenerate limits. Our studies of the thermal and 
adiabatic indices, and the speed of sound in hot and dense matter for conditions of relevance to 
core-collapse supernovae, the thermal evolution of neutron stars from their birth and mergers of 
compact binary stars reveal that substantial variations begin to occur at sub-saturation densities 
before asymptotic values are reached at supra-nuclear densities. 
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I. INTRODUCTION 

The modeling of core-collaspe supernovae, neutron 
stars from their birth to old age, and binary mergers of 
compact stars requires a detailed knowledge of the equa¬ 
tion of state (EOS) of matter at finite temperature. For 
use in large-scale computer simulations of these phenom¬ 
ena, the EOS is generally rendered in tabular forms as 
functions of the baryon density n, temperature T, and 
the electron concentration Y e = n e /n, where n e is the 
net electron density in matter. Examples of such tab¬ 
ulations can be found, e.g., in Refs. m- Entries in 
such tables include thermodynamic state variables such 
as the free energy, energy per baryon, pressure, entropy 
per baryon, specific heats, chemical potentials of the var¬ 
ious species and their derivatives with respect to num¬ 
ber densities, etc., The calculation of a thermodynam¬ 
ically consistent EOS over a wide range of densities (n 
of 10 -7 to 1-2 fm -3 ) and temperatures up to 100 MeV 
involves a detailed examination of inhomogeneous phases 
of matter (with neutron-rich nuclei, pasta-like geometri¬ 
cal configurations, leptons) at sub-nuclear densities and 
low enough temperatures, as well as homogeneous phases 
at supra-nuclear densities with possible non-nucleonic de- 
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grees of freedom (Bose condensates, strangeness-bearing 
particles, mesons, quarks, and leptons). 

Constraints on the EOS are largely restricted to zero- 
or low-temperature matter from laboratory experiments 
involving nuclei and neutron star observations. The 
various sources from which such constraints from ex¬ 
periments and astronomical observations arise are high¬ 
lighted below beginning with those at zero-temperature. 

Decades of laboratory experiments involving stable 
and radio-active nuclei have given us a wealth of data 
about nuclear masses, charge radii, neuron skin thick¬ 
nesses, nucleon effective masses, giant resonances, dipole 
polarizabilities, etc., close to the nuclear equilibrium den¬ 
sity of no — 0.16 fm -3 and isospin asymmetry a = 
(N — Z)/(N + Z) ~ (0.35 — 0.4), N and Z being the 
neutron and proton numbers of nuclei. Through mea¬ 
surements of the collective flow of matter, momentum 
and energy flow, heavy-ion collisions in the range E\ & o/A 
= 0.5-2 GeV have shed light on the EOS of matter up to 
about 3no- Astronomical observations of neutron stars 
that have central densities and isospin asymmetries sev¬ 
eral times larger than those of laboratory nuclei have 
begun to compile accurate data on neutron star masses, 
rotation periods and their time derivatives, estimates of 
radii, cooling behaviors that shed light on nucleon su¬ 
perfluidity, etc. Combined with microscopic and model 
calculations of dense nuclear matter, and the structural 
aspects and thermal evolution of neutron stars, the EOS 
of cold dense matter is beginning to be pinned down al¬ 
beit with error bands imposed by laboratory and astro¬ 
nomical data. 



2 


In specific terms, several properties of nuclei, extrap¬ 
olated to bulk matter, have yielded values for quanti¬ 
ties that characterize the key properties of isospin sym¬ 
metric nuclear matter such as the equilibrium density 
no, energy per particle Eq at no, compression modu¬ 
lus A' 0 = 9nl(d 2 E/dn 2 )\ no , and the Landau effective 
mass to*( no), with small enough one-sigma errors, al¬ 
though smaller errors would be desirable. For isospin 
asymmetric bulk matter, nuclear data have yielded val¬ 
ues for the bulk symmetry energy S v = S 2 (no) = 
(l/2)(dE(n, a)/da 2 )\ a=0 , where now the neutron excess 
parameter a = {n n — n p )/{n = n n + n p ), n n and n p being 
the neutron and proton densities in bulk matter, respec¬ 
tively. Percentage wise, the one-sigma error on S v is 
larger than those of the key parameters that characterize 
symmetric nuclear matter. Additionally, constraints on 
the derivatives of the symmetry energy at nuclear density, 

L v = 3n 0 {dS2{n)/dn)\ no and K v = 9n"o(d 2 ^(n) /dn 2 )\ no 
have also emerged, albeit with one-sigma errors that are 
somewhat large. 

Heavy-ion collisions in the energy range Ei a b/A = 
0.5 — 2 GeV have shed light on the EOS at supra-nuclear 
densities (up to 3-4 no) through studies of matter, mo¬ 
mentum, and energy flow of nucleons (for a clear exposi¬ 
tion, see [HI]). Collective flow in heavy-ion collisions has 
been characterized by (i) the mean transverse momen¬ 
tum per nucleon (p x )/A versus rapidity y/y pro j [llj . (ii) 
flow angle from asphericity analysis I12j . (iii) azimuthal 
distributions [133, and (i v ) radial flow [14]. Detection of 
neutrons being more difficult than that of protons, flow 
data gathered to date are largely for protons and for col¬ 
lisions of nuclei in which the isospin asymmetry is not 
large. Confrontation of data with theoretical calculations 
have generally been performed using Boltzmann-type ki¬ 
netic equations. One such equation for the time evolu¬ 
tion of the phase space distribution function f(r,p,t) of 
a nucleon that incorporates both the mean field U and 
a collision term with Pauli blocking of final states is (see 
Ref. [12] for a lucid account): 

% + ^pU- V r /- V r I7 • V p /= 

f- f 7 ^ 77", dtTjVJV f-Q . . 

- J d p 2 d P 2 > dll—j^-v 12 d [P + P 2 ~Pv -P 2 ') 

X [// 2 ( 1 - fv){l ~ fa) - ~ /)(1 - / 2 )] . (1) 

In general, the mean field U, felt by a nucleon due to the 
presence of all surrounding nucleons, depends on both 
the local density n and the momentum p of the nucleon. 
Operationally, U is obtained as the functional derivative 
of the energy density TL of matter at zero temperature: 
U(n,p) = SH / 5n and serves as an input on the left hand 
side of Eq. ([!]). Note that the evolution of the local 
density n is governed by the off-equilibrium evolution of 
f(r,p,t), which is driven both by spatial and momen¬ 
tum gradients of U on the one hand and by hard col¬ 
lisions on the other. The collision integral on the right 
hand side above features the relative velocity tq 2 and 
the nucleon-nucleon differential cross section dcr NN /d£l. 


Thus, Eq. ([l]) contains effects due to both soft interac¬ 
tions (the left hand side) and hard collisions (the right 
hand side) albeit at a semiclassical level insofar as phase 
space distribution functions are evolved in time clas¬ 
sically using the Boltzmann-Uehling-Uhlenbeck or the 
Vlasov-Uehling-Uhlenbeck approximation (see Ref. [10 
for original references, and a clear exposition) instead of 
the full quantum evolution of wave functions. 

Early theoretical studies that confronted data used 
isospin averaged nucleon-nucleon cross sections and mean 
fields of symmetric nuclear matter. The lesson learned 
was that much of the collective behavior observed in ex¬ 
periments stems from momentum dependent forces at 
play during the early stages of the collision nansm 
The conclusion that emerged from several studies was 
that as long as momentum dependent forces are employed 
in models that analyze the data, a symmetric matter 
compression modulus of ~ 240±20 MeV, as suggested by 
the analysis of the giant monopole resonance data EDI 
[22] , fits the heavy-ion data as well [23] ■ 

The prospects of rare-isotope accelerators (RIA’s) that 
can collide highly neutron-rich nuclei has spurred further 
work to study a system of neutrons and protons at high 
neutron excess [21M2B1 . Generalizing Eq. ([!]) to a mix¬ 
ture, the kinetic equation for neutrons is 

^ + V p A • V r /„ - V r U ■ Vp/„ = J n = J ni , (2) 

i=n,p 

where J n describes collisions of a neutron with all other 
neutrons and protons. A similar equation can be written 
down for protons with appropriate modifications. On 
the left hand side of each coupled equation the mean 
field U = U(n n , n p ;p ) depends explicitly on the neutron- 
proton asymmetry. The connection to the symmetry en¬ 
ergy arises from the fact that U is now obtained from 
a functional differentiation of the Hamiltonian density 
1-L(n n ,n p ) of isospin asymmetric matter. Examples of 
such mean fields may be found in Refs. [23H27]. Isospin 
asymmetry can be expected to influence observables such 
as the neutron-proton differential flow and the ratio of 
neutron to proton multiplicity as a function of transverse 
momentum at mid-rapidity. Experimental investigations 
of these signatures await the development of RIA’s at 
GeV energies. 

On the astrophysical front, precisely measured neu¬ 
tron star masses and radii severely constrain the EOS of 
isospin asymmetric matter. The recently well-measured 
neutron star masses 1.97 ± 0.04 M 0 [28| and 2.01 ± 
0.04 M 0 [25] have served to eliminate many, but not all, 
EOS’s in which a substantial softening of the neutron- 
matter EOS occurred due to the presence of Bose con¬ 
densates, hyperons, quark matter, etc., 30]. A precise 
measurement of the radius of the same neutron star for 
which a mass has been well measured is yet lacking, but 
reasonable estimates have been made by analysis of X- 
ray emission from isolated neutron stars, intermittently 
quiescent neutron stars undergoing accretion from a com- 
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panion star, and neutron stars that display type I X-ray 
bursts from their surfaces (3TU33| . 

A common feature shared by all of the EOS’s in Refs. 
mm is that at zero-temperature they would fail to re¬ 
produce heavy-ion flow data. The reason is in both non- 
relativistic and relativistic mean field-theoretical mod¬ 
els, the momentum dependence of the mean-field would 
lead to a linear dependence on the energy as shown 

in Refs. [13 [EMIl and [33 133 and hence at odds 

with optical model fits to nucleon-nucleus scattering data 
[13 H3 Eg 133- Potential models with finite-range inter¬ 
actions (e.g., @5j) solve this problem B3, but with at¬ 
tendant changes in their thermal properties compared to 
zero-range models as we will show in this work. A similar 
resolution in the case of field-theoretical models requires 
further work and is not addressed here. 

One of the objectives in this paper is to inquire whether 
the EOS’s that have successfully explained heavy-ion flow 
data are able to support the largest well-measured neu¬ 
tron star mass and also if the radii of 1.4 M 0 stars are 
in accord with bounds established by analyses of cur¬ 
rently available X-ray data. In order to address this issue, 
we examine models in which finite range interactions are 
employed as they yield a momentum dependence of the 
mean field U that differs significantly at high momenta 
from that of zero-range Skyrme-like models. Specifically, 
we use the MDI model of Das et al. [35] in which the 
momentum dependence was generated through the use of 
a Yukawa type finite-range interaction to match the re¬ 
sults of microscopic calculations and to reproduce optical 
model fits to nucleon-nucleus scattering data. A revised 
parametrization of the MDI model (labeled MDI(A)) was 
found necessary to support a 2 solar mass neutron star. 
We note that Danielewicz et al, show the EOS for pure 
neutron matter in Fig. 5 of Ref. [231 . but the maximum 
mass of a beta-stable neutron star was not quoted. 

Another objective of the present work is to exam¬ 
ine the thermal properties of models with finite range 
interactions and contrast them with those of a zero- 
range Skyrme model. For this purpose, we have chosen 
the SkO' parametrization m of the zero-range Skyrme 
model that yields nearly identical zero-temperature prop¬ 
erties such as the energy per particle, pressure, etc., as 
the MDI(A) model. However, the single particle po¬ 
tentials U(n n ,n p ,p) differ substantially between these 
models, the MDI(A) model being constant at high mo¬ 
menta in contrast to the quadratic rise of the SkO' model. 
Consequently, the neutron and proton effective masses, 
to* and to*, also exhibit distinctly different density de¬ 
pendences although the isospin splitting of the effective 
masses is similar in that to* > to* in neutron-rich matter 
in accordance with Brueckner-Hartree-Fock (BHF) and 
relativistic BHF calculations m- Insofar as the effective 
masses control the thermal properties, attendant differ¬ 
ences in several of the state variables become evident. 
Several of the Skyrme parametrizations exhibit a rever¬ 
sal in behavior of the isospin splitting so that to* < to* 
in neutron-rich matter (see, for example, Ref. {423 for 


a compilation, and [53] )• This has led us to establish 
conditions on the strength and range parameters of the 
Skyrme and MDI models in which to* > to*. The nu¬ 
merical results of the models studied here are supple¬ 
mented with analytical results in the limiting cases of 
degenerate and non-degenerate matter both as a check of 
our numerical evaluations and to gain physical insights. 
For use in astrophysical simulations, a detailed analy¬ 
sis of hot and dense matter containing leptons (electrons 
and positrons, and in some cases trapped neutrinos) and 
photons is also performed. Inhomogeneous phases be¬ 
low sub-nuclear densities and exotic phases of matter at 
supra-nuclear densities are not considered here as they 
fall beyond the scope of this work, but will be taken up 
in a subsequent study. 

The organization of this paper is as follows. Section |TT| 
provides an overview of the role of thermal effects in core¬ 
collapse supernovae, evolution of proto-neutron stars and 
binary mergers of compact stars. In Sec. |III[ we describe 
the finite-range (and hence momentum dependent) and 
zero-range models used to study effects of finite temper¬ 
ature. For the extraction of effects purely thermal in ori¬ 
gin, the formalism to evaluate the zero temperature state 
variables (energies, pressures, chemical potentials, etc.) 
of both these models is also presented in this section. 
Section IV presents an analysis of the zero-temperature 
results based on different parametrizations of the models 
chosen. Thereafter, differences in the momentum depen¬ 
dence of the single particle potentials between these two 
models are highlighted. Particular emphasis is placed on 
the density and isospin dependences of the Landau effec¬ 
tive masses which mainly control the thermal effects to 
be discussed in subsequent sections. In Sec. [Vj results of 
the exact, albeit numerical, calculations of the thermal 
effects are presented. This section also contains compar¬ 
isons with analytical results in the limiting cases of degen¬ 
erate and non degenerate matter. Thermal and adiabatic 
indices, and the speed of sound in hot and dense matter 
of relevance to hydrodynamical simulations of astrophys¬ 
ical phenomena involving supernovae and compact stars 
are discussed in Sec. [VH Our conclusions are in Sec. 
|VII[ Formulas that are helpful in computing the various 
state variables of the MDI model at zero temperature are 
collected in Appendix [A] In Appendix [Bj the analytical 
method by which the non degenerate limit of the MDI 
model is addressed is presented. Details concerning the 
evaluations of the specific heats at constant volume and 
pressure for the MDI and Skyrme models are given in 
Appendix [Cj 


II. THERMAL EFFECTS IN ASTROPHYSICAL 
SIMULATIONS 

The effects of temperature in astrophysical simulations 
involving dense matter are most visible in gravitational 
collapse supernovae, the evolution of proto-neutron stars, 
and mergers involving neutron stars - either neutron star- 
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neutron star (NS-NS) or black hole-neutron star (BH-NS) 
mergers. 

A. Thermal effects in supernovae and 
proto-neutron stars 

The effects of temperature in supernovae and in proto¬ 
neutron stars remain largely unexplored in detail, but 
due to the fact that for the most part matter in such 
environments is degenerate, uncertainties in the thermal 
aspects of the equation of state (EOS) do not play major 
roles in the early core-collapse phase. For example, max¬ 
imum central densities at bounce are practically indepen¬ 
dent of the assumed EOS (see, e.g., 0 )- The evolution 
of proto-neutron stars formed following core-collapse will 
be more sensitive to thermal effects, as temperatures be¬ 
yond 50 MeV are reached in the stellar cores and specific 
entropies of order lOfce are reached in the stellar man¬ 
tles |44l 115] . The maximum proto-neutron star mass and 
the evolution towards black hole formation will be depen¬ 
dent upon thermal effects. While the relative stiffness of 
the EOS (defined through the incompressibility at satu¬ 
ration) largely controls the timescale for black hole for¬ 
mation, thermal effects are important in determining the 
highest central densities reached after bounce. Reference 
[8] found that the thermal behavior was more important 
than incompressibility in this regard, largely because of 
thermal pressure support in the hot, shocked mantles of 
proto-neutron stars. In situations in which black holes do 
not form and a successful explosion ensues, binding en¬ 
ergy is largely lost due to neutrino emission and heating 
in the early evolution of proto-neutron stars, due to neu¬ 
trino downscattering from electrons, and the total neu¬ 
trino energy dwarfs the thermal energy reservoir. Nev¬ 
ertheless, the temperature at the neutrinosphere, which 
is located in semi-degenerate regions, will be sensitive to 
thermal properties of the EOS. It remains largely un¬ 
explored how the resulting neutrino spectrum, including 
average energies and emission timescales, depend on ther¬ 
mal aspects of the EOS. 


B. Thermal effects in mergers of binary stars 

Thermal effects are not expected to play a major role 
in the evolution of inspiralling compact objects up to 
the point of merger. However, the evolution of the post¬ 
merger remnant and some of the mass ejected could be 
significantly affected by thermal properties of matter. 
Perhaps the most significant recent development is the 
emergence of a standard paradigm concerning mergers of 
neutron stars. This has been triggered by the discovery of 
pulsars with approximately two solar masses, and strong 
indications that even larger mass neutron stars exist from 
a series of studies of the so-called black widow and red- 
back pulsar systems (for a review, see [33] I. Most neutron 
stars in close binaries have measured gravitational masses 


in the range 1.3 to 1.5 M 0 . The gravitational mass of the 
merger remnant will be less than twice as large, due to 
binding and due to the ejection of mass. It is unlikely 
that mass ejection will amount to more than a few hun¬ 
dredths of a solar mass, but binding energies will absorb 
larger masses. The binding energy fraction, the relative 
difference between baryon and gravitational masses, can 
be expressed by a relatively universal relation, i.e., in¬ 
dependent of the neutron star equation of state (EOS), 
involving only mass and radius [46] : 

BE N B 

- = (0.60 ±0.0 5)r A_, (3) 

where /3 = GM/(Rc 2 ). EOSs capable of supporting 
2M 0 maximum masses have the property that, for inter¬ 
mediate mass stars, the radii are nearly independent of 
the mass. Furthermore, a concordance of experimental 
nuclear physics data, theoretical neutron matter stud¬ 
ies, and astrophysical observations suggests this radius 
is about R = 12 ± 0.5 km [55] . Thus, two equal-mass 
stars with gravitational masses of 1.3M 0 will have a to¬ 
tal baryon mass of 2.87M 0 . In a merger event with no 
mass loss, a remnant of gravitational mass 2.28M 0 would 
form assuming its radius is also 12 km, representing an 
additional mass defect of O.O5M 0 relative to the initial 
stars. Should its radius decrease to about 10 km, its grav¬ 
itational mass would be approximately 1.9M 0 , and the 
additional mass defect would steepen to O.43M 0 . This 
gravitational mass could well be below the cold maximum 
mass. Repeating the above estimates for two 1.5M 0 grav¬ 
itational mass stars, we find a combined gravitational 
mass in the range 2.16M 0 -2.59M 0 , with mass defects 
larger than for the previous case. Therefore, it seems 
likely that the merged object will be close to its cold 
maximum mass. 

Studies show that the merged star will be rapidly ro¬ 
tating, and the rotation may be highly differential [47 V 
l49j . Uniform rotation can increase the maximum mass 
by a few tenths of percent (see, e.g., Ref. [50]), and 
differentially rotating objects can support further mass 
increases. In all likelihood, the merged object will be 
metastable. This possibility is enhanced if the stellar 
core is surrounded by a nearly Keplerian disc, a config¬ 
uration with an even larger metastable mass limit. If 
the merged remnant mass is above the cold gravitational 
maximum mass, but less than what rotation is capable 
of supporting, it is said to be a “hypermassive” neutron 
star, or HMNS. For a uniformly rotating star, the max¬ 
imum equatorial radius increase due to rapid rotation is 
about 50%, with a smalller change in the polar radius. 
The average density of an HMNS in that case would be 
less than half of its non-rotating value. As a result, ther¬ 
mal effects can be expected to play a much larger role in 
the stability of an HMNS than for ordinary neutron stars 
in which higher degeneracies exist. The maximum mass 
of the HMNS will decrease with time due to loss of ther¬ 
mal energy from neutrino emission and from loss of an¬ 
gular momentum due to uniformization of the differential 
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rotation mm- Early calculations, for example those of 
S3 EH], showed that collapse of the metastable HMNS 
to a black hole was induced by dissipation of differen¬ 
tial rotation and subsequent loss of angular momentum. 
In contrast, {HJ] argued that thermal effects are much 
more important than rotation in determining the stabil¬ 
ity and eventual collapse of an HMNS. In either case, the 
HMNS lifetime will crucially depend on the relative dif¬ 
ference between the HMNS mass and the value for the 
cold maximum mass for the same number of baryons. 
The HMNS lifetime, which can range from 10 ms to sev¬ 
eral seconds [3S] is potentially measurable through the 
duration of short gamma-ray bursts associated with neu¬ 
tron star mergers or the duration of gravitational wave 
signals from these mergers. Such a measurement there¬ 
fore has the potential of illuminating the EOS. 

Kaplan et al. [53] have studied the thermal enhance¬ 
ment of the pressure comparing two tabulated hot EOSs, 
LS220 [55 and Hshen [53), often used in merger simu¬ 
lations. Generally, thermal enhancements above 3 times 
the nuclear saturation density p s ~ 3 x 10 14 g cm -3 are 
less than about 5%, but at p s /2>-p s /2 the pressure is 3 
(for Hshen) or 5 (for LS220) times the cold value. As we 
will see, this difference is related to the behavior of the 
nucleon effective mass. In the stellar envelope, for densi¬ 
ties from 10 12 to 10 14 g cm -3 , the thermal enhancements 
are even larger, ranging from a factor of 10 (Hshen) to 20 
(LS220) for a thermal profile in which the average tem¬ 
perature in this density range is about 10 MeV, similar 
to that found in the merger simulations from )49j . who 
employed the Hshen EOS. In the envelope, temperatures 
are high enough such that nuclear dissociation is virtually 
complete, and the thermal pressure differences can once 
again be traced to the behavior of the nucleon effective 
mass. 

Thermal pressure alone is capable of increasing the 
maximum mass more than 10% above the cold, catalyzed 
value for a given EOS [56] . For constant temperature 
T = 50 MeV profiles studied by [53] , the maximum mass 
of a hot Hshen star was increased by 15%, while that of a 
hot LS220 star was increased by approximately 5%. For 
thermal profiles similar to the simulations of [49) , the in¬ 
creases were more modest: 3.5% and 1.8%, respectively. 

Hot rotating configurations can show either an in¬ 
crease or decrease in mass limits relative to cold rotat¬ 
ing configurations [53). When the central densities are 
less than about 10 15 g cm -3 , thermal effects increase the 
masses supported at the mass-shedding rotational limit, 
and this effect reverses at higher densities. However, the 
rotational frequency and maximum mass at the mass- 
shedding limit always decrease due to thermal support: 
the mass-shedding limit is very sensitive to the equatorial 
radius, which increases with thermal support. Thermal 
effects essentially disappear in determining the maximum 
masses of extremely differentially rotating configurations 
because the outer regions are largely Keplerian and there¬ 
fore centrifugally supported. 

The analysis by [533 concludes that thermal pressure 


support plays little role during the bulk of the evolu¬ 
tion of an HMNS, but contributes to the increase of its 
lifetime by affecting its initial conditions. As hot con¬ 
figurations with central densities less than about 10 15 g 
cm -3 support larger masses than cold ones with the same 
central density, those remnants with lower thermal pres¬ 
sure need to evolve to higher central densities to achieve 
metastability. 

In the case of BH-NS mergers, although the remnant 
will always involve a black hole, simulations indicate the 
temporary existence of a remnant disc. Although the 
disc is likely differentially rotating and largely supported 
by centrifugal effects, and its evolution controlled by an¬ 
gular momentum dissipation, thermal effects will be im¬ 
portant for dissipation due to various neutrino processes 
and emission mechanisms. 

Many early simulations of BHNS mergers employed 
cold F-law EOSs with thermal contributions that were, 
unfortunately, thermodynamically inconsistent. For ex¬ 
ample, it has sometimes been assumed that 

P = *n r + 3k ^ n +f(T) aT \ e = P/(T - 1), (4) 

where k and T are constants and f(T ) is a temperature- 
dependent factor reflecting the fraction of relativisitic 
particles in the gas, ranging from 1 at low temperatures 
to 8 at high temperatures. However, not only should / 
also be density-dependent, but also it can be shown from 
the above expression for e that, even if it is treated as 
being density-independent, 

P = nn r +aT j f{T)T 2 dT . (5) 

Obviously these two expressions for the pressure are in¬ 
compatible. Therefore, interpreting the thermal behav¬ 
iors of merger calculations with T-law EOSs is problem¬ 
atic. 

Direct comparisons of BHNS merger simulations with 
tabulated, temperature-dependent EOSs have yet to be 
made. Individual simulations with tabulated EOSs have 
been performed in [53 with the Hshen EOS and in [58) 
with the LS220 EOS. Most focus has been on the proper¬ 
ties of the remnant disc formed from the disrupted neu¬ 
tron star which survives on timescales ranging from tens 
of ms to several s. Typically, densities in the early evolu¬ 
tion of the disc range from l-4x 10 11 g cm -3 , proton frac¬ 
tions are around Y e ~ 0.1, and specific entropies range 
from 7-9 ks- As is the case with HMNS evolutions, the 
thermal properties control neutrino emissions and the ul¬ 
timate disc lifetimes. 

III. MODELS WITH FINITE- AND ZERO- 
RANGE INTERACTIONS 

A. Finite range interactions 

We adopt the model of Das et al. [23) who have gen¬ 
eralized the earlier model of Welke et al. ca to the case 
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of isospin asymmetric nuclear matter. In this model, ex¬ 
change contributions arising from finite range Yukawa 
interactions between nucleons give rise to a momentum 
dependent mean field. BUU calculations performed with 
such a mean held have been able to account for data 
from nuclear reactions induced by neutron rich nuclei 
[55] , Recently, Ref. m has reported results of thermal 
properties of asymmetric nuclear matter. The model’s 
predictions for the structural properties of neutron stars, 
and thermal effects for conditions of relevance to astro- 
physical situations have not been investigated so far and 
are undertaken here. Explicitly, the MDI Hamiltonian 
density is given by [23] 

R = ^(T n + t p ) + V(n n ,n p ,T), (6) 


where 


n * = J d3pi and Ti= J d3 Pi 


are the number densities and kinetic energy densities of 
nucleon species i = n,p, respectively. The potential en¬ 
ergy density V = V(n n ,n p ,T ) is expressed as 
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Above, no — 0.16 fm ~ 3 is the equilibrium density of 
isospin symmetric matter. We will discuss the choice of 
the strength parameters Ai, A 2 , B, y , Ci , C u , the pa¬ 
rameter (T that captures the density dependence of higher 
than two-body interactions, and of the finite range pa¬ 
rameter A in subsequent sections. For simplicity, the 
finite range parameter A is taken to be the same for both 
like and unlike pairs of nucleons, but the strength param¬ 
eters, Ci and C u , are allowed to be different. Accounting 
for 2 spin degrees of freedom, the quantities 


fi{n,Pi) 


2 

( 2irh ) 3 


fpi ~ 


2 1 
( 2nh ) 3 1 + 


(9) 


are the phase-space distributions of nucleons in a heat 
bath of temperature T, having an energy spectrum 


e 


Pi 


2 &h &h 

dri dni 


2m 


Y Ui{n n , n p , pi) 


( 10 ) 


where 


Ui{Ui, Uj,Pi) = —{Ui Y Hj) + —{Hi - Uj) 
no no 


B 


ni + rij 

n 0 


ni + n, 


1 Y 


i -y 

2 


(7—1 

(7+1 


ni + nj 


(7 — 1 V ni — nj 


Y Ri(ni,n.j,pi) 

— Ui{ni,nj) + Ri(rii, nj,pi) , i 4 j (H) 


is the single-particle potential (or the mean field). Above, 
Ui represents contributions arising from the densities 
alone. The momentum dependence is contained in 


R. l (n l ,n j ,p l ) 


2 Q 2 f 3 , _/+ 

n 0 (27 rh) 3 J 1 (jh- 


2 C u 2 

no (27 rh) 3 


d 3 Pj 



(.12) 


At finite temperature, the determination of Ri(rii,nj,Pi) 
requires the knowledge of R, t {rii, n 3 ,p)) for all values of 
p'. As in Hartree-Fock theory, a self-consistency condi¬ 
tion must be fulfilled; this is achieved through an iterative 
procedure as in Ref. m- The initial guess is supplied 
by the zero-temperature Ri(ni, Uj,Pi) (analytical expres¬ 
sions are given in the next section) which may be used in 
Eq. (H) to obtain a starting chemical potential and 


energy spectrum ef\pi). These are then used in Eqs. 
(10) and ([7]) to obtain Ri{ni,nj,Pi) from Eq. (12). This 
in turn leads to e^(pi) which is used in Eqs. (10) and 
Q to find Upon repetition of the cycle, convergence 
is achieved in five or less iterations for most cases. The 
ensuing chemical potentials /ii, for a given density, tem¬ 
perature and composition, are then used in the standard 
statistical mechanics expression for the entropy density: 



[f Pi 1 nf Pi Y (1 - f Pi ) ln(l - f Pi )\. 


(13) 


By integrating this expression twice by parts, s for the 
MDI model takes the form 



Y 2 


d 3 Pi , 
( 2 t tH) 3 Ip ‘ 


Ri{Pi ) Y 


Pi dRj(pi) 

3 d Pi 


(14) 


The pressure is acquired through the thermodynamic 
identity 


P = — £ + Ts + PiHi. 


(15) 
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where the energy density e = R. The result in Eq. (14) 
enables pressure to be cast in the form m 

P = (n n + n p ) 2 + E- (n n - n p ) 2 

ZTLq Z77-0 


aB ( n n + n P Y +1 


a + 1 


1 - 2 / 


d 3 pi 


(2ttH) 3 p% 1 3 m 


fr> , i E 


Vi dRi Ri i, 


The same expression can also be obtained from Eq. (3.6) 
of Ref. [IT]. 

We follow Ref. [61 j to express the specific heat at con¬ 
stant volume as 


C v = 


1 de 
n dT 


(17) 


and the specific heat at constant pressure as m 

r (^1 ) 2 

1 \dT\n) 


C P = Cv + 


OP I 
dn It 


(18) 


By performing a Jacobi transformation to the variables 
p, and T, 


de 

Cv= m ; 


C P = Cv + 


de 


dn 1 

dfi 

T 

dT 1 fj. 


dn 



dfi 

T 


/ 

ap | 

dn | \ 

f 3PI 

dn It 

3t| m \ 

9T ,! 

dn 

| J 

V p 

dpi 

It / 


dp | 
dn It 

O" I 
dn It 


(19) 


( 20 ) 


respectively. 

The calculation of Cy and Cp for the MDI model in¬ 
volves some intricacies not encountered in the zero-range 
Skyrme-like models as U(rii,n p ,p) in Eq. (10) depends 


only on the nucleon densities and the temperature. Con¬ 
sequently, derivatives in the above equations must be 
evaluated with some care as described in Appendix [C| 


B. Zero temperature properties 

The MDI Hamiltonian density can be written as the 
sum of terms arising from kinetic sources, R &, density- 
dependent interactions, Rd, and momentum-dependent 
interactions, R m : 


( 21 ) 


At T = 0, 


_ 1 1 1 

k 2 to Tn + Tp 2m 5 tt 2 H 3 

Rd = T-n 2 + ^n 2 (l - 2x) 2 

zno 277-o 


(.PFn 


Pf p ) 


B n a+l 


R „ 


a + 1 n\ 

= C i (I 

n 0 


o 


[l-y(l-2x) 2 } 


2Cuj 

J VP) i ± np 

n 0 


with 


PFi = (37t 2 mh 3 ) 1/3 

8 tt 2 A 2 

(2irh) 

4A 

e^-{PFi — PFj 
4A 


( 22 ) 

(23) 

(24) 

(25) 


hj = ^PFiPFjiPpi + PFj) - 


PFiPFjR 2 

3 


(PFi ~ PFj) arctan 


-y (Pt* + Pt 7 ) arctan 


PFi - PFj 
PFi + PFj 

A 


x In 


A 4 ; {p 2 Fi +P 2 Fj)R 2 _ ( P 2 Fi -PFj? 
12 + 2 4 

(PFi + PFj) 2 + A 2 


(PFi - PFj) 2 + A 2 _ 


(26) 


We note that Eq. (3.5) in Ref. (35], which agrees with 
Eq. (26 1 above for isospin symmetric matter and for pure 


neutron matter, must be corrected to properly account 
for properties of bulk matter with intermediate isospin 
content. 

The energy per particle, the pressure, and the chemical 
potentials are obtained from the relations 


E=* 

n 


P = n^-R. 
on 


and pn = — (27) 

uni 


where n = n n + n p is the total baryon number density. 

Symmetric nuclear matter properties at the saturation 
density no = 0.16 frn -3 such as the compression modulus 
A'o, the symmetry energy S v as well as its slope L v and 
curvature K v are obtained from 


Kq = I\(n = no, x = 1 / 2 ) = 9ip 


d 2 R 


dn 2 


(28) 


n=riQ , x— 1/2 


S v = S 2 (no ) ; S 2 = 7T- 
L v = L 2 (n 0 ) = 3 n 0 


1 d 2 R 


8 n dx 2 
dS 2 


®=l /2 


K v = K 2 {n 0 ) = 9n 0 


dn 

2 d 2 S 2 
dn 2 


(29) 

(30) 

(31) 


The proton fraction x is defined as x = n p /(n n + n p ). 
The T = 0 single-particle energy spectrum is 


p ? 

£Pi — "h Ui(m, nji hi) 
2 m 


'H = 7^/e + Hd + . 


(32) 









































with 


form, the Skyrme Hamiltonian density reads as 


Ui (jlli , Tlj , ki ) — Ai 


fa + rij) + fa - rij) 


+ B 


n 0 

Hi + rij 
n 0 


,, 1 , 1 

ft = --+ --T, 


?r 0 


i -y 


( 7—1 


2?77.. 

+ n(r n + T p ) 


2m v p 


rii — rij 
Hi + rij 


1 + 


(7+1 

2 / n, + rij 


cr — 1 \ rii — n 


•3 / J 


+ ( 


'Tri'T^n i TpTlp 


fafafafa) 

ti A , 


T Bii{nii ki) T Rij fa, ki) , 2 7^ J 
n 0 n 0 

(33) 


W 1 
4 fa 
(1 


£2 


zi 


4 fa 

10 fa) "'fa (fa) « + + 

11 (‘ + ^) " 2 - il (I+^) 


where 


(38) 


fa fa,fa 


A 3 f 2+p 


Explicit forms of the single-particle potentials for the 
SkO' Hamiltonian (with 1) = rii/n, i = n,p ) are 


4f 2 H 3 1 A 
- 2 


A 


(A 2 +fa -fa) 

+ -;^fa-hr 


arctan I Pl +PFj ) - arctan < Pi PFj 


2 A pi 


A 

fa +pfa 2 + A 2 ' 
.fa ~PFj) 2 + A 2 _ 


[fan, fc) = (AA + YjX 2 )nk 2 + {X 1 + X 2 )n + A^- 
+ 2n(A3 + YiX 4) + {(2 + e)A5 

+ [2Y i + e(Y i 2 + Y j 2 )]X 6 }; faj, 


(39) 


(34) where 


from which the nucleon effective masses can be derived: 
* t ddn 


ITli =p Fi 


Explicitly, 


nii = p F i 


dpi 


-1 


PFi 


1 + 


m dUj 

PFi dpt 


(35) 


Xl = 4L 
to 


h (l + y) +^(l + f) 

*2 ( ^ + X 2 J - t\ ( ~ + Xi 




PFi 


A, = fa 1 + fa 


12 


fafa 


*o fa , 

= ~2 2 +*° 


A 6 = —— f- + x 3 
6 12 \ 2 3 


(40) 


/ PFi ! 2 Ci 

d Ru 

! 2 C u 

dRij 

\ 

_i masses are 


l rn tiq 

dpt 

PFi ^ 

dpi 

PFi , 

) ’ < = 

/HI 

1 + yy (fa + YiX 2 )n 






(36) m 


From Eq. (35), the density-dependent Landau effective 

(41) 


-1 


dRi 


dpi 


A 2 


PFj 


2ii 2 H 3 pf 


1 - 


(A 2 +p 2 F1 +P 2 Fi) 


x In 


4 PFiPFj 

C PFi + PFj) 2 + A 2 ' 

(PFi - PFj) 2 + A 2 _ 


For completeness, explicit expressions for 
E, P , /j,; , A, fa, L 2 , A 2 and m* as functions of 
the nucleon densities and isospin content are collected 
in Appendix A. 


The derivations of the various state variables, nu¬ 
clear saturation properties and thermal response func¬ 
tions proceed as previously described for MDI. For details 
, .of evaluating the thermal state variables for Skyrme-like 
models, we refer the reader to a recent compilation of 
formulas and numerical methods in Ref. [63| . For nu¬ 
merical values of the various strength parameters above, 
we choose the SkO' parametrization of Ref. [4(3] which 
will be given in a subsequent section. 


C. Zero-range Skyrme interactions 

For comparison with the results of the finite range 
model discussed above, we also consider the often stud¬ 
ied zero-range model due to Skyrme [62] . In its standard 


IV. RESULTS FOR ZERO TEMPERATURE 

In this section, we consider the zero temperature prop¬ 
erties of the finite-range and zero-range models discussed 
in the previous section. We begin with the MDI(O) 
parametrization of the model of Das et al. (24] so that 
its characteristics extended to neutron-star matter, not 
considered previously, may be assessed. 
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A. MDI models for isospin-asymmetric matter 

Table [I] lists the various parameters employed in the 
MDI(O) model of Ref. [H]. In Table [Tl] , we list the 
characteristic properties of this model at the equilibrium 
density of isospin symmetric nuclear matter. Also in¬ 
cluded in this table are values of the various physical 
quantities (the last three rows) accessible to laboratory 
experiments for small isospin asymmetry. Note the fairly 
good agreement with experimental determinations of the 
various quantities. 


Parameter 


Value 

Parameter 


Value 

Ai - 

T08.28 

MeV 

y 

0 


A 2 

-12.30 

MeV 

Ci 

-11.70 

MeV 

B 

106.35 

MeV 

C u - 

T03.40 

MeV 

a 

4/3 


A 

263.04 

MeV 


TABLE I. Parameter values for the MDI(O) Hamiltonian den¬ 
sity of Das et al. [24,. The dimensions are such that the 
Hamiltonian density is in MeV fm -3 . 


Property Value Experiment Reference 


[MDI(O)] 


77,0 (fm d ) 

0.16 

0.17±0.02 

[MH53 

E 0 (MeV) 

-16.10 

—16±1 

E3E2 

A'o (MeV) 

212.4 

230±30 

[23123 



240±20 

m 

mo/m 

0.67 

0.8T0.1 

[69jEQ! 

S v (MeV) 

30.54 

30-35 

EDEZ1 

L v (MeV) 

60.24 

40-70 

[ID EH 

K v (MeV) 

-81.67 

— 100T200 

m 


TABLE II. Entries in this table are at the equilibrium density 
no of symmetric nuclear matter for the MDI(O) model [24] . 
Eo is the energy per particle, A'o is the compression modulus, 
mj/m is the ratio of the Landau effective mass to mass in 
vacuum, S v is the nuclear symmetry energy parameter, L v 
and K v , are related to the first and second derivatives of the 
symmetry energy, respectively. 


In Table III the structural properties of neutron stars 
built using the EOS of charge neutral and beta equili¬ 
brated matter from the MDI(O) model are summarized. 
The predicted maximum mass falls slightly short of the 
largest well-measured mass. The radii of the maximum 
mass and 1.4 M 0 stars are in reasonable agreement with 
their current estimations from X-ray data. A notewor¬ 
thy feature is the central baryon chemical potential of 
the maximum mass star, fi c ~ 1.9 GeV, which is below 
the limit of 2.1 GeV set by the maximally compact EOS 
of a neutron star derived in Ref. [3D]. Note that the 
density n a , at which the EOS violates causality, that is, 
the squared speed of sound c 2 = dP/de > 1 (a common 
feature of potential models), lies well above the central 
density of the maximum star as is the case for all models 
discussed in subsequent sections. 


Property 

Value 

[MDI(0)] 

Observation 

Reference 

M max(Mo) 

1.884 

2.01±0.04 

m 

Rmax (km) 

9.84 

11.Oil.0 

m 

n c (fm -3 ) 

1.3065 



e c (MeV fm -3 ) 

1703.3 



P c (MeV fm" 3 ) 

761.22 



Ate (MeV) 

1886.3 



n a (fm -3 ) 

1.883 



Pi.4(km) 

11.77 

11.5T0.7 

m 

n c (fm -3 ) 

0.5782 



e c (MeV fm" 3 ) 

599.19 



P c (MeV fm" 3 ) 

91.784 



Ate (MeV) 

1195.0 




TABLE III. M m ax is the maximum neutron star mass pre¬ 
dicted by the MDI(O) model and Rmax the radius associated 
with it. Other entries are the central density n c , energy den¬ 
sity e c , pressure P c , and the chemical potential p c for both 
the maximum mass and i.4 M@ configurations, and the den¬ 
sity n a at which the EOS becomes acausal. The value for the 
radius of a 1.4 M@ neutron star is given by Pi.4. 


B. Revised parameterization of the MDI model 


In this section, we provide a revised set of param¬ 
eters (see Table IV ) for the MDI model so that the 
isospin symmetric and asymmetric properties at the nu¬ 
clear matter equilibrium density are closer to the experi¬ 
mentally derived mean values than given by the MDI(O) 
parametrization and the neutron star maximum mass 
also comes close to the recently observed 2 solar mass. 
We note that this model also allows us to constrain the 
single particle potential U(n,p) to match optical model 
fits to data. We have used U(no,p = 0) = —74.6 MeV 
and U(no,p 2 /(2m ) = 303.1 MeV) = 0 in determining the 
constants of the model as the variational Monte Carlo cal¬ 
culations of Wiringa in Ref. m suggest. The asymptot¬ 
ically flat behavior of {/(?ro,p) — 30.6 MeV with p arises 
naturally from effects of the finite range interaction in 
this model. Before presenting the results of the MDI(A) 
model, we consider a zero-range model for purposes of 
comparison. 


Parameter 

Value 

Parameter Value 

Ai 

-69.4758 MeV 

y -0.0327929 

A 2 

-29.2241 MeV 

Ci -23.0576 MeV 

B 

100.084 MeV 

C u -105.885 MeV 

a 

1.36227 

A 420.864 MeV 


TABLE IV. Parameter values for the MDI(A) Hamiltonian 
density. The dimensions are such that the Hamiltonian den¬ 
sity is in MeV fm -3 . 
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C. Parameters of a prototype Skyrme model 


As an example of a zero-range model, we have cho¬ 
sen to work with the Skyrme model of Ref. [40] known 
as SkO', the parameters of which are listed in Table 
[Vl For an apposite comparison, the parameters of the 
MDI(A) model were tuned so that its energy per particle 
vs baryon density closely matches that of the SkO' model 
for both symmetric matter and pure neutron matter. 



U 

Xi e 

0 

-2099.419 MeV fm 3 

-0.029503 1/4 

1 

301.531 MeV fm 5 

-1.325732 

2 

154.781 MeV fm 5 

-2.323439 

3 

13526.464 MeV fm 3(1+<0 

-0.147404 


TABLE V. Parameter values for the Skyrme Hamiltonian den¬ 
sity SkO' used in this work. The dimensions are such that the 
Hamiltonian density is in MeV fm -3 . 


Comparison of MDI(A) and SkO' models 


Attributes of the MDI(A) and SKO' models at their re¬ 
spective equilibrium densities of symmetric nuclear mat¬ 
ter are presented in Table m The resulting structural 
aspects of neutrons stars from these two models are pre¬ 
sented in Table |VII| These results indicate the nearly 
identical nature of the two models at zero temperature. 
As with the MDI(O) model, the central baryon chemical 
potentials for these models also lie below the value for the 
maximally compact EOS. A closer examination of the in¬ 
nards of these two models is taken up in the subsequent 
sections. 


Property 

Value 

[MDI(A)] 

Value 

[SkO'] 

Experiment 

Reference 

no (fm -3 ) 

0.160 

0.160 

0.17T0.02 

EHSZ] 

E 0 (MeV) 

-16.00 

-15.75 

—16±1 

[6UE3 

K 0 (MeV) 

232.0 

222.3 

230±30 

[SEMI] 




240±20 

[61 

mj$/m 

0.67 

0.90 

0 .8T0.1 

El EH 

S v (MeV) 

30.0 

31.9 

30-35 

eh El 

L v (MeV) 

65.0 

68.9 

40-70 

E0ES] 

K v (MeV) 

-72.0 

-78.8 

-100T200 

ESI 


TABLE VI. Same as Table [!!] but for the MDI(A) and SkO' 
models. 


We wish to add that neutron star maximum masses 
in excess of 2 M 0 can also be obtained from a 
reparametrization of the MDI model, but at the ex¬ 
pense of losing close similarity with the results of the 
SkO' model. An illustration is provided with K () = 260 
MeV, L v = 70 MeV, and K v = —50 MeV, while keeping 
the other saturation properties the same as for MDI(A). 
This was achieved through the choice of the constants 


Property 

Value 

[MDI(A)] 

Value 

[SkO'] 

Observation Reference 

M max(Mo) 

1.9725 

1.9600 

2.01T0.04 

m 

Rmax (km) 

10.20 

10.13 

11 .Oil.0 

m 

n c (fm -3 ) 

1.2065 

1.2233 



e c (MeV fm -3 ) 

1573.4 

1595.7 



P c (MeV fm -3 ) 

718.90 

739.42 



Ate (MeV) 

1900.0 

1908.9 



n a (fm -3 ) 

1.64 

1.678 



Pi.4(km) 

12.21 

12.17 

11.5T0.7 

“Esr 

n c (fm -3 ) 

0.5126 

0.5234 



e c (MeV fm -3 ) 

526.36 

533.06 



P c (MeV fm -3 ) 

75.802 

78.421 



ATc(MeV) 

1174.7 

1168.3 




TABLE VII. 
models. 


Same as 


III 


but for the MDI(A) and SkO' 


A i = -39.0752, A 2 = -27.4916, B = 69.6838, a = 
1.56497, y = -0.105539, Cj = -26.95, C u = -101.993, 
and A = 420.864, the units of these constants being the 
same as in Table [T] Results for the structural proper¬ 
ties of a neutron star in beta-equilibrium are: M max = 
2.153 M 0 , Rmax = 10.58 km, and R\a = 12.13 km. At 
the edges of the 1-cr errors of the empirical saturation 
properties at the nuclear equilibrium density, it is not 
difficult to raise the maximum mass well above 2 Mg. 


D. Single particle potentials 


In this section, we present results of the single particle 
potentials for the MDI(A) and SkO' models from Eqs. 
<[33|) and ([39| , respectively, and contrast them with those 


from the microscopic calculations of Refs. [74j and m- 
Figure shows the neutron single particle potentials 
for MDI(A) and SkO' as functions of momentum for se¬ 
lect baryon densities at zero temperature. Results shown 
are for pure neutron matter (x = 0) [Fig. 1(a)] and for 
isospin asymmetric matter with x = 0.2 [Fig. 1(b)]. No¬ 
tice that the results for MDI(A) tend to saturate at large 
momenta for both proton fractions owing to the loga¬ 
rithmic structure of Eq. (341. The SkO' model, how¬ 


ever, in common with most Skyrme models, exhibits a 
quadratic rise with momentum. This latter feature is 
also present in the results of the Akrnal, Pandharipande 
and Ravenhall (APR) model [76] in which the hamilto- 
nian density of the many-body calculations of Akrnal and 
Pandharipande m is parametrized in Skyrme-like fash¬ 
ion. For both MDI(A) and SkO' models the effect of a 
finite proton fraction [Fig. 1(b)] is more pronounced at 
low momenta for which the single particle potential be¬ 
comes more attractive relative to that for pure neutron 
matter. 

In Fig. [2] the neutron single particle potentials vs 
momentum for the MDI(A) model from Eqs. (33) 
(solid curves) and for the SkO' model from Eq. (34) 


(dashed curves) are compared with the variational Monte 
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FIG. 1. Neutron single particle potentials vs momentu m at zero temperature for the baryon densities marked. Results for the 
MDI(A) model (solid curves) are from from Eqs. (101, (331, and (34|. Those for the SkO 7 model (dashed curves) are from Eq. 
(139}. (a) Results are for proton fraction x = 0. (b) Same as (a) but for x = 0.2. 


Carlo results (solid curves marked with asterisks) of Ref. 
m using the UV14-TNI interaction and the Bruekner- 
Hartree-Fock results (dash-dotted curves) of Ref. [75] 
with the inclusion of three-body interactions (labeled 
BHF-TBF). Results are for symmetric nuclear matter at 
about one, two and three times the nuclear matter equi¬ 
librium density in figures (a), (b), and (c), respectively. A 
curve-to-curve quantitative comparison between the re¬ 
sults of models used in this work and those of UV14-TNI 
and BHF-TBF models is not appropriate because the 
saturation properties of the latter models differ signifi¬ 
cantly from those of the former ones. Specifically, for the 
UV 14-TNI model, E 0 = -16.6 MeV at n 0 = 0.157 fm -3 , 
with K 0 = 260 MeV, whereas for the BHF-TBF model, 
Eq = —15.08 MeV at rig = 0.198 fm -3 , with Kq = 207 
MeV. In this work, the single particle potential was de¬ 
signed to mimic closely the behavior of the UV14-TNI 
model. However, the qualitative trends - power-law rise 
vs logarithmic rise - at high momenta in other micro¬ 
scopic models are worth noting as discussed below. 

At the densities shown, and at high momenta, there 
is good agreement between the results of the MDI(A) 
model and those of the UV14-TNI model by design. As 
mentioned in the introduction, this saturating behavior 
at high momenta is in accord with the analysis of optical 
model fits to nucleon-nucleus scattering data. The other 
two models, SkO' and BHF-TBF, rise quadratically with 
momenta. In the case of the SkO' model, this behavior 
ensues from the zero-range approximation made for nu¬ 
clear interactions. In the case of the BHF-TBF calcula¬ 
tions, the quadratic rise with momentum stems from the 
similar behavior of U (n, k ) chosen for purposes of conver¬ 
gence during numerical calculations. Such a behavior of 
the single particle potential, even with the symmetric nu¬ 


clear matter compression moduli around 240 MeV, leads 
to nucleon collective flows that are larger than those ob¬ 
served in heavy-ion collisions na ushui- We note that 
Danielewicz reaches similar conclusions with a different 
parametrization U (n, p ), which also saturates at high mo¬ 
menta (see. Fig. 17 of Ref. [T8]). 


300 

250 

200 
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100 

50 

0 

-50 

-100 


300 
250 
200 
150 
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50 
0 

-50 
-100 

1234 1234 12345 
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FIG. 2. Comparison of neutron single particle potentials vs 
momentum from Eqs. ([33} and ([34} with the variational 
Monte Carlo results of Ref. (74] using the UV14-TNI inter¬ 
action and Bruekner-Hartree-Fock results of Ref. ES3 with 
the inclusion of three-body interactions (labeled BHF-TBF). 
Results shown are for zero temperature symmetric nuclear 
matter at the indicated baryon densities. 
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FIG. 3. Landau effective masses of the neutron (solid curves) and proton (dotted curves) scaled with the vacuum nucleon mass 
vs baryon density n for the marked values of the proton fraction x. (a) Results are for the MDI(A) model from Eqs. (361-( 371. 
(b) Same as (a), but for the SkO' model from Eq. (351. 


E. Isospin dependence of effective masses 

The single particle potentials discussed above facilitate 
the calculation of the Landau effective masses of nucleons 


according to Eq. (35). For the MDI(A) and SkO' models, 


explicit expressions as functions of density and proton 
fraction were given in Eqs. ( |37| ) and (41), respectively. 

In Fig. [3j the neutron and proton effective masses 
scaled with the vacuum nucleon mass are shown as a 
function of baryon density for select proton fractions for 
the MDI(A) and SkO' models in panels (a) and (b), re¬ 
spectively. Although these two models yield similar prop¬ 
erties for most observables for symmetric nuclear matter 
at the equilibrium density n 0 , the effective masses are sig¬ 
nificantly different - mj/m = 0.67 (0.9) for the MDI(A) 
(SkO') model (see Table VI). The density dependence of 
m n,p/ m also differs significantly between the two mod¬ 
els - a logarithmic decrease in the MDI(A) model vs a 
[1 + (constant) • ?r] _1 decrease with density in the SkO' 
model. Effects of isospin content as it varies from that of 
symmetric nuclear matter (x = 0.5) toward pure neutron 
matter (x —> 0) are qualitatively similar, but quantita¬ 
tively different with MDI(A) producing a significantly 
larger change compared with SkO'. It is worthwhile to 
note, however, that several parametrizations of Skyrme 
interactions exist in the literature which yield a larger 
variation of to* p /m with varying x than is present in the 
SkO' model, although the logarithmic decline with den¬ 
sity of the MDI models would be absent in all of them. 

A noteworthy feature of the results in Figs. [3](a) and 
(b) is that m* > to* for all densities as x moves from its 
symmetric matter value of 0.5 to 0, the value for pure 
neutron matter. The cause for this behavior may be 
traced to the values of strength and range parameters 


that govern the behavior of effective masses with pro¬ 
ton fraction. For example, if we require the condition 
to* > to* to be satisfied for the MDI models, Eq. (37) 
implies that to leading order in n and for all x 


Cl - C u > 0, 


(42) 


independent of the finite range parameter A. Next to 
leading order in n and for all x, the condition becomes 


Ci 


-Cu 


1-f ttf+4 


(*) C^) 


( n \ 

( mC u \ 

\ n °J 

V A 2 ) 


>0, (43) 


where pf above is the ferrni momentum of symmetric nu¬ 
clear matter. As (20/3 )(p_f/A) 2 > 4(py/A) 2 , the condi¬ 
tion becomes a Ci > C u with a < 1. These conditions are 
met for the MDI(0) (and for all the MDI models in Ref. 
1211) and MDI(A) models which ensures that to* > to* 
for all x in the range 0.5 — 0. 

For Skyrme interactions, Eq. (41) implies that to* > 


m* in neutron-rich matter as long as 

ti(1 + 2xi) > t 2 (l + x 2 ), 
which is satisfied by the parameters in Table V. 


(44) 


When the conditions in Eqs. (42), (43), and (44), are 
not satisfied, a reversal in the behavior of neutron and 
proton effective masses occurs, that is, to* < to* as x 
moves away from 0.5 toward 0. Notwithstanding this 
behavior, the properties of symmetric nuclear matter, 
symmetry energy attributes, and collective excitations 
of nuclei have been well described. Additionally, the re¬ 
quirement that the EOS of neutron-star matter is able 
to support stars of 2M 0 has also been met. To illus¬ 
trate this, we show in Table VIII parameters of an MDI 
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model and those of a Skyrme model (SLy4 [75]) in Table 
IX whose predictions at the nuclear equilibrium density 
are summarized in Table X. The corresponding results 
for neutron star structure are given in Table XI. The 
Landau effective masses for the neutron and proton that 
exhibit this reversal in behavior as a function of x are 
shown in Figs. |4j)a) and (b). Although we have not 
done an exhaustive search, we have verified that several 
parametrizations of Skyrme interactions, GSkI, GSkll, 
KDEOvl >79, '80] to name a few, exhibit this reversal in 
behavior of neutron and proton effective masses as a func¬ 
tion of x. 


Parameter 

Value 

Parameter 

Value 

Ai 

290.942 MeV 

y 

1.47111 

A 2 

709.676 MeV 

a - 

-525.764 MeV 

B 

151.01 MeV 

c u 

-0.0587396 MeV 

a 

1.23671 

A 

1315.20 MeV 


TABLE VIII. Parameter values for the MDI(B) Hamiltonian 
density. The dimensions are such that the Hamiltonian den¬ 
sity is in MeV fm -3 . 


Property 

Value 

[MDI(B)] 

Value 

[Sly4] 

Observation Reference 

M max(Mo) 

2.0295 

2.037 

2.01T0.04 

m 

Ftmax (km) 

10.027 

9.797 

11.Oil.0 

m 

n c (fm -3 ) 

1.1768 

1.2062 



e c (MeV fm -3 ) 

1539.5 

1598.8 



P c (MeV fm" 3 ) 

749.91 

856.45 



Ate (MeV) 

1945.4 

2035.54 



Pi.4(km) 

11.41 

11.12 

11.5T0.7 


n c (fm -3 ) 

0.5110 

0.5413 



e c (MeV fm -3 ) 

522.41 

555.99 



P c (MeV fm" 3 ) 

77.390 

86.391 



Ate (MeV) 

1173.8 

1186.604 




TABLE XI. Same as[m] but for the MDI(B) and Sly4 models. 


els with similar splitting exhibit quantitative differences 
only. The isospin splitting of the effective masses in the 
MDI(A) model is in agreement with the above micro¬ 
scopic models. 


i 

u 

Xi 

e 

0 

-2488.91 MeV fm 3 

0.834 

1/6 

1 

486.82 MeV fm 5 

-0.344 


2 

-546.39 MeV fm 5 

-1.000 


3 

13777.0 MeV fm 3(1 + e ) 

1.354 



TABLE IX. Parameter values for the Skyrme Hamiltonian 
density Sly4 EB3- The dimensions are such that the Hamilto¬ 
nian density is in MeV fm -3 . 


Property 

Value 

[MDI(B)] 

Value 

[Sly4] 

Experiment 

Reference 

no (fm -3 ) 

0.160 

0.160 

0.17T0.02 

SMHSZ] 

E 0 (MeV) 

-16.00 

-15.97 

—16±1 

[Ml ill 

K 0 (MeV) 

239.0 

229.9 

230T30 

240T20 

EU [ 22 ] 
[MS 

m.o/m 

0.67 

0.69 

0.8T0.1 

ESI [70] 

S v (MeV) 

30.0 

32.0 

30-35 

eh [72] 

L v (MeV) 

52.0 

46.0 

40-70 

EH ED 

K v (MeV) 

-81.0 

-119.7 

-100T200 

m 


F. Energy, pressure and chemical potentials 


The energy per baryon E = 1~L/n from the MDI(A) and 
SkO' models is presented as functions of baryon density 
and proton fraction in Fig. [5] Explicit expressions for 7~L 
are in Eq. (241 for the MDI model and in Eq. (381 for 


the SkO' model. For all proton fractions ranging from 
pure neutron matter to symmetric nuclear matter shown 
in Fig. [5ja) there is little difference between the energies 
of the two models, the energy for the SkO' model being 
slightly larger than that of MDI for pure neutron matter 
at all densities. The inset in Fig. [5])b) shows the energy 
per baryon of nucleons vs baryon density for charge neu¬ 
tral matter in beta equilibrium. The two models yield 
nearly the same energy at most densities with only small 
differences around the nuclear matter density. 

In terms of the neutron excess parameter a = ( n n — 
n p )/n = 1 — 2x , the PNM to SNM energy difference can 
be expressed as 


AE = E(n, a) — E(n, 0) = Si(n)c 


(45) 


1 = 2 , 4 ,... 


TABLE X. Entries in this table are at the equilibrium den¬ 
sity no of symmetric nuclear matter for the MDI(B) and Sly4 
models. Eq is the energy per particle, Kq is the compression 
modulus, mo/m. is the ratio of the Landau effective mass to 
mass in vacuum, S v is the nuclear symmetry energy parame¬ 
ter, L v and K v , are related to the first and second derivatives 
of the symmetry energy, respectively. 

It is worthwhile noting that in the microscopic 
Brueckner-Hartree-Fock and Dirac-Brueckner-Hartree- 
Fock calculations that include 3-body interactions, to* > 
to* in neutron-rich matter |81I484| . Their density de¬ 
pendences, while akin to those of Skyrme-like mod- 


where 

; l = 2,4,... (46) 

cu—0 

In Fig. |6])a) , AE of the MDI(A) model is compared 
with its symmetry energy S -2 (n) which provides a reason¬ 
able approximation. Contributions of the higher order 
terms stemming from kinetic and momentum dependent 
interaction terms (denoted by Sik and Si m , respectively) 
up to l = 6 are shown Fig. §b). The convergence of 
AE — Si to zero is shown in Fig. [6](c). At the nu¬ 
clear equilibrium density, S 2 falls short of AE by a little 


1 dE(n,a) 
1 l\ da 1 
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E 
* E 



n (fm -3 ) 


n (fm -3 ) 


FIG. 4. Landau effective masses of the neutron (solid) and proton (dotted) scaled with the vacuum nucleon mass vs baryon 
density n for the marked values of the proton fraction x. (a) Results are for the MDI(B) model from Eqs. (36 1-(|37|. (b) Same 


as (a), but for the Sly4 model from Eq. (351. Note the reversal in behavior of the effective masses when compared to the results 
in Figs. [3]) a) and (b). 



n (fm' 3 ) 


FIG. 5. (a) Zero temperature energy per particle E = T-L/n 
vs baryon number density n at the p roton fra ctions x shown 
for the MDI(A) model from Eqs. (22 1 -( 241 and the SkO 7 
model from Eq. ( |38| ). (b) E vs n for matter with varying 
proton fractions xp determined from charge neutrality and 
beta-equilibrium. 



n (fm' 3 ) 


FIG. 6. (a) The difference A E = E(n, x = 0) — E(n, x = 1/2) 
vs n and the symmetry energy S 2 vs n. (b) Higher order 
contributions of Si to A E. (c) Convergence of A E — Si up 
to 6th order terms in Si. 


over an MeV. The inclusion of higher order terms in Si 
drastically reduces the difference. Figure [7] shows simi¬ 
lar results for the SkO' model. Relative to the MDI(A) 
model, the largest qualitative difference in seen in the 
density dependence of Sp m . However, the convergence of 
A E — Si to zero is similar for both models. 

In Fig. §a), we show the pressures exerted by nucleons 
resulting from the MDI(A) and SkO' models as functions 
of density for symmetric nuclear and pure neutron mat¬ 
ter. The results in Fig. 8jb) correspond to charge neu¬ 
tral matter in beta equilibrium. Both models have nearly 


identical pressures regardless of the proton fraction with 
small differences at the highest densities shown. 


The density dependence of the neutron and proton 
chemical potentials from the MDI(A) and SkO' mod¬ 
els are presented in Figs. and (b), respectively. 

The M DI( A ) che mical potential was calculated using 
Eqs. (A7)-(A10). For all proton fractions presented 
(x = 0.1, 0.3 and 0.5), the neutron chemical potentials 
of the two models agree at all densities. Small differ¬ 
ences between the two models are observed in the proton 
chemical potentials with the largest difference occuring 
at low proton fractions and large baryon densities. Fig- 
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0.5 

3 
2 
1 
0 

0 0.3 0.6 0.9 0.2 0.4 0.6 0.8 

n (frrf 3 ) 

FIG. 7. Same as Fig. [ 6 ] but for the SkO' model. 




FIG. 8 . (a) Zero temperature pressure P vs baryon density 
n at t he p r oton fractions shown for the MDI(A) model from 
Eqs. (All-(A 6 l and the SkO 7 model from Eq. (271. (b) P 
vs n for matter with varying proton fractions xp determined 
from charge neutrality and beta-equilibrium. 


ure |9jc) shows the difference between the neutron and 
proton chemical potentials (p = p n — p p ) vs density for 
the MDI(A) and SkO' models, respectively. As there is 
little difference between the neutron chemical potentials 
of the two models, the small differences between the p's 
stem from differences in the proton chemical potentials. 
The differences in p occur mainly for low proton fractions 
and high densities. 


G. Properties of cold-catalyzed neutron stars 


try energies in the two models (characterized by similar 
values of L at no). 

The central densities of stars with different masses are 
of relevance in the long-term cooling of neutron stars j85| . 
For up to a million years after their birth, neutron stars 
cool primarily through neutrino emission. As pointed 
out in Ref. [85], the density dependence of more 

precisely A E(n), determines the threshold densities for 
the rapid direct Urea processes n O p + e~ + D e and 
n •<-> p+p~ + D p to occur. In charge neutral (n e - + n p - = 
n p ) and beta-equilibrated (p = p e - = p ^-) neutron- 
star matter, the threshold densities for these direct Urea 
processes are determined when the triangular inequalities 

PFi + PFj > PF k , (47) 

where i,j and k are p, e~ (ft - ), and n , and cyclic per¬ 
mutations of them, are satisfied (neutrinos at the rel¬ 
evant temperatures contributing little to the momen¬ 
tum balance conditions). For the direct Urea process 
involving electrons (and muons), the threshold density 
is 0.67(0.93) fm -3 for the MDI(A) model. Thus, stars 
with masses M > 1.7Mq will undergo rapid cooling from 
neutrino emission from the process involving electrons 
while those with M > 1.9Mq will receive an equal and 
additional contribution to neutrino emissivity from the 
process with muons. The corresponding numbers for the 
SkO' model are 0.59(0.92) fnU 3 and 1.5(1.9)M 0 , respec¬ 
tively. 


V. THERMAL EFFECTS 

Finite-temperature effects can be studied by isolating 
the thermal part of the various functions of interest de¬ 
fined as the difference between the finite-T and T = 0 
expressions for a given quantity X: 

X th = X(rii, nj,T) - X(ni,rij, 0 ). (48) 

For Skyrme-like models, we refer the reader to Sec. V of 
Ref. [63] where explicit expressions and numerical notes 
are provided for evaluations of the thermal state variables 
and response functions. Thus, our discussion below will 
focus on evaluations of the thermal state variables for the 
MDI models only. 

Suppressing the symbols denoting dependences on the 
baryon density n and proton fraction x for brevity, the 
thermal energy is given by 


In Figs. [M a ) and (b), we show the gravitational 
mass vs radius ( Mg vs R ) and mass vs central den¬ 
sity (Mq vs n c ) of cold-catalyzed neutron stars for the 
MDI(A) and SkO' models. Not unexpectedly, the re¬ 
sults are very nearly the same as the P vs n relations 
in neutron-star matter are nearly identical for the two 
models [see Fig. ] 8 ](b)]. The close to vertical rise of mass 
with radius is due to the similar behaviors of the symme- 


E th = E(T) - E( 0) 


(49) 


where E(T) = H(T)/n is calculated using Eqs. @-([ 12 ]) 
in Sec. II-A, and E( 0) = T-L(0)/n is obtained from Eqs. 
[24])-[26]) in Sec. II-B. 


Likewise, the thermal pressure is 


Pth = P(T) - P( 0), 


(50) 
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FIG. 9. [(a) an d (b )1 The neutron and proton chemical potentials vs baryon density n for the MDI(A) [Eqs. (AT I- (A101] and 
the Sk(y [Eq. (271] models for values of x shown, (c) ft = p n — vs n. 



2 

1.5 

1 

0.5 

0 


FIG. 10. (a) Mass vs radius and (b) mass vs central den¬ 
sity of cold-catalyzed neutron stars for the MDI(A) and SkO' 
models. 


where P(T) and P( 0) are calculated using Eq. (16 1 in 
Sec. II-A, and in Appendix | A 1| respectively. Details re¬ 
garding the numerical evaluation of the chemical poten¬ 
tials at finite T are described in Sec II-A in connection 
with Eq. (12 1 and the discussion thereafter. Explicit ex¬ 
pressions for evaluating the chemical potentials at T = 0 
are collected in Appendix |A 2| Utilizing these results, the 
thermal chemical potentials are given by 


particle spectrum in Eq. (10 1 can be written as 


2 2 

£p = !ro + U ( n,P ^ = irn + R ( n ’P'’ T ) + u i n ) > ( 52 ) 


where the explicit p-dependence arising from finite-range 
interactions is contained in R(n,p;T) and U(n) is a 
density-dependent, but p-inclependent term from inter¬ 
actions that are local in space. For a given n at fixed 
T, the chemical potential p is determined from Eq. ([T]) 
using an iterative procedure as outlined earlier. At finite 
T, the term R acquires a T-dependence (see, e.g., Fig. 


11 1 which has been explicitly indicated. We can write 


£p ~ M = 2m + R ^ n ' p ’ T '> ~ \P~ u ( n )\ 

= +R(n,p;T) - R(n,p = 0 ;T) 

-[fi-{U(n) + R(n,p = 0-,T)}], (53) 


where in the last step we have isolated the p-independent, 
but n- and T-dependent part i?.(n, p = 0; T) and grouped 
it with U(n). Note that i?(n,p = 0; T) provides an n- and 
T- dependent pedestal for the momentum dependence of 
R(n,p;T). As a result, 


€p T ^ ^ { 7^ + R (n, P', T ) — R{n,p = 0; T)} -77 

(54) 

with 


Pth = p{T) - p(0). 


(51) 


The entropy per baryon S = s/n issues from Eq. (JT3J) of 
Sec. Il-a, the specific heats at constant volume and pres¬ 
sure, Cy and Cp, from Eqs. and ( |I8| ), respectively, 
of the same section. 

The question arises as to what plays the role of the 
degeneracy parameter for a general momentum depen¬ 
dent interaction. The exposition below is for a 1- 
component system. Generalization to the 2-component 
case is straightforward. For the MDI model, the single 


p= ^[p-{U{n) + R(n,p = 0-T)}}, (55) 

which serves as the degeneracy parameter. The term 
R(n,p;T) when combined withp 2 /(2m) generates an ef¬ 
fective mass in*. (For free gases, R{n,p\T = 0) = 0 and 
U{ri) = 0 so that 77 = p/T.) 

For Skyrme interactions, R(n,p;T) = /3 x np 2 /(2m), 
where /3 is an interaction strength-dependent constant. 
Combined with the p 2 /(2m) term, the p-dependent part 
of e p can be written as p 2 / (2m*) with m* = m( 1 + pn) -1 






















17 


being a density-dependent, 
mass. Thus, 


but T-independent, effective 


(see Fig. it is useful to study the T = 0 values for 
orientation. In this case (with f p / = 1 and y = p' / A), 


p = —— ^ ^ for Skyrme interactions . (56) 

For a general R(n 7 p), more complicated dependencies 
of m* with n ensue including a T-dependence as for the 
MDI models with Yukawa-type interactions. 

Utilizing 

/j,(n,T = 0) = E F = ^ + R(n,p F -,T = 0) +U(n) (57) 
2m 

(as follows from Hugenholtz-Van Hove’s theorem) we also 
have 


V= f {p(n,T = 0) + p th (n 7 T) - \U(n) + R(n,p = 0;T)]} 
= f lff + [R{n,p F ;T = 0)-R(n 7 p = 0- 7 T] + p th (n,T) 


where the thermal contribution to the chemical potential, 
Pth, is as defined in Eq. (511. Above, p F is the Fermi 


momentum at T = 0 determined using p F = fi,(37r 2 n) 1//3 . 
Recall that the number of particles is fixed be it at T = 0 
or T 0. Consequently, p F serves as a useful reference 
momentum with respect to which effects of finite T (for 
which states with p < p F and p > p F up to oo are all 
occupied with finite probabilities) can be gauged. For 
example, 




■ P n) + p t h(n,T) 


1 

T 


Pt 

2m* 


+ p t h(n,T) 


(59) 


for Skyrme interactions. For MDI models, an equivalent 
expression can be derived but it features additional den¬ 
sity dependencies and agrees with Eq. (59) only when 


the p-dependence is truncated at the quadratic term. 

The sign of p depends on the magnitudes of the various 
terms in Eq. (58). The first two terms are positive, and 
generally, p t h < 0 for all degeneracies (see Fig. 12). 


In the degenerate limit, the sum of the positive terms 
dominates over p, t h rendering p positive, whereas in the 
non-degenerate limit the magnitude of pth exceeds the 
summed positive terms making p negative. 

The contribution of R(n 7 p = 0 ;T) to Eq. (55) is sig¬ 


nificant as an analysis for the MDI model used here il¬ 
lustrates. For like particles, 


on 

R(n,p = 0-,T) = — 2 
n 0 


d 3 p' 


(2nh) 3 fp ' 1 + (£)2 ’ 


1 


(60) 


where C = Ci for pure neutron matter (PNM) and C = 
Ci + C u for symmetric nuclear matter (SNM). Although 
temperature corrections are significant at low momenta 


R(n,p = 0;T = 0) = — (j) 4 
n 0 \nJ tt z 


dy 


1 +y 2 


( 61) 

where yf = p F / A. The integral above is easily performed 
with the result 


R(n 7 p = 0; T = 0) = 


2C 


n 0 


4- tan” 1 4 

L A A 


= 6C ( — 

n 0 


_ £ (VF_y , £ (P£_Y _ 


m u( 


(62) 


The results in Fig. |TT] are in agreement with the above 
result in the degenerate limit. Corrections from temper¬ 
ature effects would be of 0{(T/E F ) 2 ). 

An analytical result for the non-degenerate limit can 
e obtained from Eq. (B12): 


C(n\ A 2 

n,p = 0;T) = — — - 

9 \noJ Por 


1 , A 2 
- m-rw- 
p A 2 


(P + POR ) 


21 


(p-Por) 2 


(63) 

where g = 2(4) for PNM (SNM) and por = {2mT) 1 / 2 . 
From a Taylor expansion of the term in brackets, 


p —>0 


R(n,p = 0; T) = 


c ( 

r n \ 

A 2 

4poR 

9 ' 

K n 0 ) 

Por 

A2 +Por 

4 C 

/ n N 

\ 

A 2 

9 

V^Oy 

1 A 2 

+ 2mT 


(64) 


The non-degenerate results in Fig. 0 are in accord with 
expectations from the above result. 

In what follows, we compare the various state variables 
and the specific heats of the MDI(A) and SkO' models. 
Results of these calculations are presented in Figs. [TT| 
through [20] for pure neutron matter and symmetric nu¬ 
clear matter at T = 20 and 50 MeV, respectively. 


A. Results of numerical calculations 

For given n 7 x , and T, the chemical potentials p n and 
g p for the MDI(A) model can be calculated using the 
iterative procedure described earlier in connection with 
Eq. (12). This entails a self-consistent determination of 


11 


i? ra (p)and R p {p ), results for which are shown in Fig. 
for pure neutron matter and isospin symmetric nuclear 
matter at various densities, and temperatures of T = 20 
and 50 MeV, respectively. In few iterations, convergence 
is reached from the starting guess R^ at T = 0 to the 
final result RC) shown in panels (a) through (d) of this 
figure. As expected, substantial corrections to the zero- 
temperature result at low to intermediate momenta are 
evident as the temperature increases. 

The thermal chemical potentials of the neutron p n jh 
for the two models are shown in Fig. |T2] Both mod¬ 
els predict that the p n ,th ’s of symmetric matter and pure 
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FIG. 11. Momentum dependent part of the neutron single particle potential for the baryon densities, temperatures, and proton 
fractions indicated in each panel. The zeroth approximation R <0> is the zero temperature R„ gi ven by the momentum dependent 
terms in Eq. (331, whereas the final result R^' 1 is found via an iterative procedure using Eq. (121. 


neutron matter cross with the latter having a larger value 
at low densities. The density at which the crossing oc¬ 
curs is smaller in the case of MDI(A) than that of SkO'. 
Larger temperatures move the crossing to larger densi¬ 
ties. For pure neutron matter the two models predict 
nearly the same result at all densities and temperatures. 
The saturating behavior in both models is a consequence 
of progressively increasing degeneracy of the fermions 
and its onset occurs at higher densities for the higher 
temperature. 



n (fm" 3 ) 


50 
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-50 


-100 
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FIG. 12. Thermal chemical potentials vs n at the indicated 
proton fractions x and temperatures T for the MDI(A) and 
SkO' models from Eq. |t]| . 


Contours of the degeneracy parameter are shown in 
Fig. [13] While the qualitative trends are similar for the 
MDI(A) and SkO' models, quantitative differences exist 
in the degenerate regime for r) > 1. The origin of these 
differences can be traced to the different behaviors of the 


effective masses in the two models. 

The thermal energy per baryon from the MDI(A) and 
SkO' models are shown as functions of baryon density in 
Figs. 14’a) and (b). At T = 20 MeV, the MDI(A) model 


has somewhat lower values than those of the SkO' model 
beyond nuclear densities for symmetric nuclear matter 
(x = 0.5). However, results of the two models agree to 
well beyond the nuclear density for pure neutron matter 
(x = 0), significant differences occurring only for densi¬ 
ties beyond those shown in the figure. An opposite trend 
is observed at T = 50 MeV for which the two models dif¬ 
fer slightly around nuclear densities for x = 0.5, whereas 
they yield similar results for x = 0 at subnuclear densi¬ 
ties. We attribute these behaviors to the significantly dif¬ 
ferent behaviors of the effective masses (both their mag¬ 
nitudes and density dependences) in these two models 
(see Fig. [3]) as our analysis in the subsequent section, 
where analytical results in the limiting cases of degen¬ 
erate and non degenerate matter are compared with the 
exact results, shows. In order to further appreciate the 
extent to which thermal energies differ quantitatively be¬ 
tween different models, it is instructive to compare these 
results with those in Fig. 10 of Ref. |63j where pre¬ 
dictions for the APR and the Skyrme-Ska models were 
recently reported. 

The thermal free energy per baryon is shown in Fig. 15 
as a function of baryon density for the MDI(A) and SkO' 
models. Results for the two models are indistinguish¬ 
able at low densities (n < 0.01 fm -3 ). This low den¬ 
sity agreement between the two models improves with 
increasing temperature and with lower proton fractions. 
For n > 0.01 fm -3 , quantitative differences between the 
two models are due to the different trends with density 
of nucleon effective masses. These differences become in¬ 
creasingly small for all proton fractions as the limit of 
extreme degeneracy is approached at very high densities. 
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FIG. 13. Contours of the degener acy parameter r/ in the n-T plane. Results for the MDI(A) model are from Eq. (55 I and those 
for the SkO' model are from Eq. (561. 
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FIG. 14. Thermal energy per particle vs n at the indicated 
proton fractions x and temperat ure s T for the MDI(A) and 
SkO' models from Eqs. |6| and (381. 


FIG. 15. Thermal free energy for the MDI(A) and SkO' mod¬ 
els at the indicated tem pera tures an d proton fractions. Re¬ 
sults are from Eqs. (Jg]) , (381 and (131. 


In Fig. 16 we present the thermal pressures vs density. 
For both temperatures and proton fractions shown, the 
two models display similar traits in that at around nu¬ 
clear and subnuclear densities they predict similar values, 
but begin to differ substantially at supra-nuclear densi¬ 
ties. With increasing density, the thermal pressure of the 
MDI(A) model is smaller than that of the SkO' model 
chiefly due to its smaller effective mass and relatively 
flat variation with density. 


Figure [TT] shows the entropy per baryon for the two 
models. The two models agree at low densities with the 
best agreement occurring for pure neutron matter up to 
about twice the nuclear density. At larger densities the 
MDI(A) model predicts that the entropy of symmetric 
matter converges to that of pure neutron matter. This 
feature is also present in SkO' but occurs at larger den¬ 


sities than shown here. 


Isentropic contours in the n-T plane are shown in Fig. 


18 for a proton fraction of 0.1 and entropies in the range 


0.25-3. Both models show similar trends in that all con¬ 
tours rise quickly until around no/2, beyond which only 
a moderate increase in the temperature is observed. For 
each entropy contour, the temperature is systematically 
larger for the SkO' model when compared with that of 
the MDI(A) model. For densities larger than 2no and for 
values of entropy exceeding 1.5, the temperatures pre¬ 
dicted by both models are well in excess of 50 MeV. 


The specific heat at constant volume, Cy, is plotted 
as a function of baryon density in Fig. [19] for the two 
models. Noteworthy features at both temperatures are 
the peaks occurring at values in excess of 1.5 (the maxi¬ 
mum value characteristic of free fermi gases at vanishing 
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FIG. 16. Thermal pressure vs n at the indicated proton frac¬ 
tions x and temperatures T for the MDI(A) and SkCF models 
from Eq. fl5| ). 


Fig. 20 as a function of baryon density. The predominant 
feature for both models is the sharp rise in Cp for sym¬ 
metric nuclear matter at T = 20 MeV. This feature arises 
from the temperature being close to that for the onset of 
the liquid-gas phase transition for which dP/dn = 0. At 
high densities Cp resembles the behaviors seen for Cy 
and the entropy per baryon. 


Figure [2l] shows the ratio of specific heats Cp/Cy vs n 
for T = 20 and 50 MeV, respectively, for the two models. 
In (a), the large variations seen at sub-saturation densi¬ 
ties for values of x not too close to that of pure neutron 
matter are due to the proximity of an incipient liquid- 
gas phase transition. As will be discussed in the next 
section, Cp/Cy is closely related to the adiabatic index, 
Tg, which provides a measure of the stiffness of the equa¬ 
tion of state. In addition, it also determines the speed 
of adiabatic sound wave propagation in hydrodynamic 
evolution of matter. 



FIG. 17. Entropy per particle vs n at the indicated proton 
fractions x and temperatures T for the MDI(A) and SkO' 
models from Eq. (13 1 . 


density, which is also the case for Skyrme models) at fi¬ 
nite densities in the MDI(A) model. These peaks can 
be attributed to the momentum dependence built into 
the interaction which produces a temperature-dependent 
spectrum via R(n,p). This trait, shared with relativis¬ 
tic mean field models (although there, the T-dependence 
in the spectrum enters through the Dirac effective mass) 
[86] , has implications related to the hydrodynamic evolu¬ 
tion of a core-collapse supernova in that Cy controls the 
density at which the core rebounds. Cy decreases with 
increasing density and the magnitude of the decrease is 
larger at the lower temperature. As was the case for 
the entropy per baryon, the Cy 's of nuclear and neutron 
matter approach each other at large densities. 

The specific heat at constant pressure, Cp, is shown in 


B. Analytical results in limiting cases 

In the cases when degenerate (low T, high n such that 
T/Ep i 1) or non-degenerate (high T, low n such that 
T/Ep i 1) conditions are met, generally compact ana¬ 
lytical expressions can be derived. As the densities and 
temperatures encountered in the thermal evolution of su¬ 
pernovae, neutron stars, and binary mergers vary over 
wide ranges, matter could be in the degenerate, partially 
degenerate, or non-degenerate limits depending on the 
ambient conditions. A comparison of the exact, but nu¬ 
merical, results with their analytical counterparts not 
only allows for a check of the often involved numerical 
calculations, but is also helpful in identifying the density 
and temperature ranges in which matter is one or the 
other limiting case. Because of the varying concentra¬ 
tions of neutrons, protons, and leptons, one or the other 
species may lie in different regimes of degeneracy. 


1. Degenerate limit 


In this case, we can use Landau’s Fermi Liquid Theory 
(FLT) [87l [88] to advantage. The ensuing analytical ex¬ 
pressions highlight the importance of the effective mass 
in thermal effects. The explicit forms of the thermal en¬ 
ergy, thermal pressure, thermal chemical potentials and 
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FIG. 18. Isentropes in the n-T plane for the MDI(A) model in (a) and those for the SkO' model in (b). 




FIG. 19. [(a) and (b)] Specif ic h eat at constant volume Cy 
vs baryon density from Eq. (171 for the MDI(A) and SkO' 
models. 


FIG. 20. [(a) and (b)] Specific heat at constant pressure Cp 
vs baryon density from Eq. (181. 


entropy density in FLT are as follows: 


T 2 

E t h = — V' dint 
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(65) 

( 66 ) 

(67) 

( 68 ) 


where cq = ^ Cy- is the level density parameter. In this 

limit, to lowest order in temperature, Cy = Cp = s/n. 
The above relations are quite general in character requir¬ 
ing only the concentrations and effective masses (which, 



n (fm“ 3 ) 


FIG. 21. Ratio of specific heats at the indicated proton frac¬ 
tions and temperatures. 


in turn, depend on the single particle spectra) of the var- 
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ious constituents in matter. 


2. Nondegenerate limit 

Non-degenerate conditions prevail when the fugacities 
Zi = e p " i ^ T are small. Methods to calculate the state 
variables in this limit for Skyrme-like models have been 
amply discussed in the literature (see, e.g., Ref. f53] for a 
recent compilation of the relevant formulas) and will not 
be repeated here. For the MDI models in which the single 
particle spectrum receives significant contributions from 
momentum-dependent interactions, the analysis is some¬ 
what involved. We have developed a method involving 
next-to-leading order steepest descent calculations that 
provides an adequate description of the various state vari¬ 
ables (see Appendix [B]for details). The numerical results 
presented below for the non-degenerate limit are obtained 
employing the relations in Appendices [B] and [C] 


C. Numerical vs analytical results 


In this section, the exact numerical results of Sec. [V] 
are compared with those using the analytical formulas in 
the degenerate and nondegenerate limits described in the 
previous section. Throughout, results from the MDI(A) 
model are displayed in panels (a) and (b), whereas panels 
(c) and (d) contain results from the SkO' model in Figs. 
(23|-([27|) below. 


Figure [22] contains plots of the exact thermal chem¬ 
ical potential of the neutron and, its degenenerate and 
nondegenerate limits. The agreement between the non¬ 
degenerate limit and the exact result is significantly bet¬ 
ter for SkO' compared with MDI(A). This is best seen 
in the T = 50 MeV results for pure neutron matter for 
which the nondegenerate limit coincides with the exact 
result until about 2 n 0 for SkO' compared with MDI(A) 
which agrees only to about uq. Both models predict con¬ 
vergence between the degenerate and exact results begin¬ 
ning at around 2no for T = 20 MeV and around 4ng for 
T = 50 MeV. 

The exact thermal energy and its limits from the two 
models are shown in Fig. 23 The MDI(A) model has a 
thermal energy that agrees with its non-degenerate limit 
for similar densities compared to that for the SkO' model. 
For both models the agreement between the exact result 
and the non-degenerate limit is better at high tempera¬ 
tures and for pure neutron matter. Agreement between 
the degenerate limit and the exact solution occurs sooner 
(lower density) for MDI(A) than for SkO'. In both cases, 
the best agreement is for lower temperatures and for sym¬ 
metric nuclear matter. Note, however, that around nu¬ 
clear densities, matter is in the semi-degenerate limit as 
was the case for APR and Ska models in [53]. 

In Fig. [24[ we present the thermal pressure and its 
limiting cases. For both models the non-degenerate limit 
agrees with the exact result until about no at 20 MeV 


(panels (a) and (c)) and until 3?zo for T = 50 MeV (pan¬ 
els (b) and (d). The agreement is better at high temper¬ 
atures and for symmetric matter than for pure neutron 
matter. At both temperatures and proton fractions, the 
degenerate limits come closer to the exact results using 
the MDI(A) model. 

In Fig. [25| we present the entropy per baryon and, 
its degenerate and non-degenerate limits. The non¬ 
degenerate limit has the best agreement using SkO' for 
symmetric nuclear matter at T = 50 MeV, which ex¬ 
tends to about 3no- The range of densities over which 
the MDI(A) model agrees with the exact result is smaller 
than that for SkO'. For symmetric nuclear matter, for 
example, the agreement does not extend beyond about 
1-1.5 no even at 50 MeV. The degenerate limit coincides 
with the exact solution starting approximately around 
2no for the MDI(A) model for symmetric nuclear matter 
at T = 20 MeV. The SkO' model does notably worse 
with its best agreement not occurring until 3-4 no for 
pure neutron matter at T = 20 MeV. 

Some insight into the behaviors of the thermal vari¬ 
ables presented above can be gained from the asymp¬ 
totic behaviors of the single-particle potentials at high 
momenta in the two models which lead to two notewor¬ 
thy effects: (1) The earlier onset of degeneracy in the 
MDI model compared to that for SkO' is due to weaker 
binding at high densities, and (2) the MDI nucleon ef¬ 
fective masses that are nearly independent of density at 
high density (while qualitatively the isospin splitting is 
similar to that of SkO') cause the thermal state variables 
to exhibit less sensitivity to the proton fraction being 
changed (cf. FLT equations with dm*/dn —► 0). 

The specific heat at constant volume vs baryon den¬ 
sity and its limiting cases are presented in Figs. [26[a) 
and (b) for the MDI(A) model, while those for the SkO' 
model are in (c) and (d) of the same figure. The best 
agreement between the results of the exact and the de¬ 
generate limit calculations occurs at low temperatures 
and large densities. Although this is true of both mod¬ 
els, the agreement is better for the MDI(A) model as the 
degenerate limit comes far closer to the exact solution 
than the SkO' model. For MDI(A), the degenerate limit 
has better agreement with the exact solution for symmet¬ 
ric nuclear matter as opposed to SkO' which shows better 
agreement for pure neutron matter. The non-degenerate 
limit coincides with the exact solution only for densi¬ 
ties much less than the nuclear saturation density. The 
agreement between the non-degenerate limit and the ex¬ 
act result is best using the MDI(A) model for symmetric 
nuclear matter at high temperatures. The agreement be¬ 
tween the non-degenerate limit and the exact result using 
the SkO' model is slightly better for pure neutron matter 
and at high temperatures. 

In Fig. [27] we display the specific heat at constant 
pressure vs baryon density from the MDI(A) (panels (a) 
and (b)) and SkO' (panels (c) and (d)) models. The 
agreement between the non-degenerate limit and the ex¬ 
act solution is remarkably good using the MDI(A) for 
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FIG. 22. Neutron thermal chemical potentials [Eq. (JtJIvs n for the MDI(A) model in (a) and (b) and the SkO ; model in (c) 
and (d) compared with their limiting cases from Eqs. (671 and (B201. 




80 


60 


40 


20 


0 


FIG. 23. Thermal energy per particle [Eq. (16])] and limiting cases [Eqs. (65 I and (B31)] vs n for the temperatures and proton 


fractions shown. 
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FIG. 24. Thermal pressure [Eq. and limiting cases [Eqs. ( |66| ) and ( |B36| )] vs n for the temperatures and proton fractions 
shown. 


pure neutron matter at high temperatures and low den- sities. For the MDI(A) model the agreement extends out 
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FIG. 25. Entropy per particle [Eq. (13 1] vs n for the MDI(A) model in (a) and (b) and the SkO' model in (c) and (d) compared 
with their limiting cases from Eqs. (|68|) and (B35). 
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FIG. 26. Specific heat at constant volume CV for the MDI(A) model in (a) and (b) and the SkO' model in (c) and (d) compared 
with their limiting cases from Eq. (681 as S = Cv in the degenerate limit, and Eq. (C231. 


to about 0.4 fm~ 3 , whereas for the SkO' the agreement 
is up to 0.3 fm~ 3 for symmetric or pure neutron matter 
at high temperatures. The agreement between the de¬ 
generate limit and the exact result for Cp is best for the 
MDI(A) model for pure neutron matter at large densities 
and low temperatures. Using the SkO' model, the agree¬ 
ment between the degenerate limit and the exact solution 
is better for pure neutron matter at large densities and 
low temperatures. 

VI. THERMAL AND ADIABATIC INDICES 

The paradigm for neutron star mergers now seems to 
be that one begins with two 1.3 -1.4 Mq/Mq stars, which 
form a hyper-massive remnant stabilized against collapse 
by rotation, thermal and magnetic effects. It could be dif¬ 
ferentially rotating. Loss of differential rotation, and loss 
of thermal and/or magnetic support leads to an eventual 


collapse to a black hole. The timescale is very impor¬ 
tant, as it will have observable effects on gravitational 
wave and gamma-ray burst durations. Thermal effects 
at supra-nuclear density seem to have little effect, but 
rotational support means the average densities of the 
disc are near saturation density where the thermal ef¬ 
fects become substantial (see references below). Thus, 
the thermal support needs to be properly treated. 


A. Thermal index 

The inclusion of thermal effects in neutron star merger 
simulations is often treated using an effective thermal 
index T t h{n) defined as 

T th (n) = l + — . (69) 

£th 
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Shibata’s group commonly uses T t h to describe finite- 
temperature effects, favoring the value 1.8 [55]. Bauswein 
et al., [50] prefer the value 2.0; see also Janka, et al. En¬ 
in simulations by Foucart et al., [22> and Kaplan et al., 
[98] conditions such as S' = 5 at 10 14 g cm -3 and S = 10 
at 10 12 g cm -3 are reached in the ejecta. In these works, 
realistic EOS’s with consistent thermal treatments (LS 
[Ml or Shen EE among others) are used. An overview 
has been provided in the work of Cyrol |94] in which 
the behavior of of T t h vs. n has been presented using 
the tabulated results of non-relativistic and relativistic 
models. 

Our aim here is to provide a basis for understanding 
the behavior of r ( /, from elementary considerations. The 
results of our calculations are relevant only for densities 
and temperatures for which a bulk homogeneous phase 
will be present. Inhomogeneous phases present at sub- 
nuclear densities, and which induce large variations in 
T fh , have not been considered here as they lie beyond 
the scope of this work. 


In the degenerate limit, the FLT results in Eqs. (65) 
and (661 imply that 


r th (n) = i+ 2 -^ iainiQl 


i = n,p,e 


3 J2i a i n i 

where the level density parameters are 

non — relativistic nucleons 


a,- = 


2 pi 


rzp. _|_ m 2 

v » - 1 relativistic electrons . 


2 p 


and 


Qi = i - 


3 n dm* 
2 to* dn 


(70) 


(71) 


(72) 


The above equations highlight the role of the effective 
masses and their behavior with density. Note that rela¬ 
tivity endows non-interacting electrons with a density- 


dependent effective mass to* = E Fb 
Thus, in a pure electron gas, 


Qe 


1 p'k 

2 p 2 p + TO 


2 ’ 
e 



(73) 


which has the limit 1/2 for ultra-relativistic electrons 
(pf c TO e ) and 1 for non-relativistic electrons (p Fe •C 
m e ). These limits help to recover the well known results 
Fth = 4/3 in the former case and 5/3 in the latter (also 
easily obtained by inspecting the limits of Pth/zth in the 
two cases). 


Note that in the degenerate limit, T t h in Eq. (70) is 


independent of temperature. For pure neutron matter 
(PNM) and symmetric nuclear matter (SNM) without 
electrons, Eq. (70) reduces to the simple result 


, . 5 n dm t 

m = 3 -^*r 


(74) 


where the subscript “6” identifies the appropriate 
baryons (neutrons in PNM, and neutrons and protons 
in SNM or isospin asymmetric matter). 

For Skyrme models, the above relation is valid for all 
regions of degeneracy as Pth and £th can be written in 
terms of their ideal gas counterparts (calculated with 
to* (n) instead of to) as 


P th (n,T) = Ptf(n,T-,m*) ( 1—£ 
e th(n,T) = e\ d h (n,T-,m*) 


3 n dm* 


2 to* dn 
2 


Tjid 
r th 

r~.id 

^th 


(75) 


These results in conjunction with Eq. (69) lead to Eq. 


(74). The simple form of the effective masses in Skyrme 
models, to* = to( 1 + /3n) -1 , where the positive constant 
p depends mildly on the proton fraction (for the SkO' 
model, P lies in the range 0.523-0.724 as Y p varies from 
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0-0.5) allows us to obtain 


r th(n) 


5 fin 
3 + 1 + fin 


8 m* 
3 m 


(76) 


which establishes the T- independence and very mild de¬ 
pendence on the proton fraction. For Skyrme models 
therefore, T t h of nucleons increases monotonically from 
5/3 to 8/3 as the density increases. 

The analytical expressions for the effective masses and 
their derivatives with respect to density are more com¬ 
plicated for the MDI models than those for the Skyrme 
models (see Appendix A). However, they are easily im¬ 
plemented in numerical calculations. Results from such 
computations will be compared with the exact numerical 
calculations below. 

In the Maxwell-Boltzmann limit [that is, the non¬ 
degenerate limit to 0(z 1 )}, the thermal energy density 
and pressure of the MDI model are given, respectively, 
by 


£th — 


3 nT 


exp 


R(y/2mT) - R(V^mT) 


T 


+ ™R(V2rtT) - -T F n 

£ o 

Rfij2mT) 


Pth = nT 


2 T 


- 

5 


(77) 

(78) 


where R(p) is given by Eq. (B12) with pofi = V / 2mT 
and Tf = p F /{2m*). Exploiting the fact that in this 
regime the interactions are weak due to the diluteness of 
the system, we expand the exponential in Eq. (771 in a 


Taylor series about the zero of its argument which leads 
to 


£th — 


3 nT 


1 + 


4/3R(V2toT) - R(v'4mT) 


T 


- -Tfu. 

5 

(79) 

Note that in Eqs. ([78]) and ([79]) the leading terms are 
proportional to nT, the interaction terms to n 2 /T 2 (ap¬ 
proximately) and the T = 0 terms to n 5 ^ 3 . Thus in the 
limit of vanishing density, the ratio Pth/tth goes to 2/3; 
consequently, T f h approaches 5/3 as expected for a non- 
relativistic gas. 

In order to appreciate the role of electrons in the be¬ 
havior of T th v s n, we first show in Fig. [28] results with 
nucleons only for three different models. Proton frac¬ 
tions and temperatures are as noted in the figure. Re¬ 
sults for the non-relativistic models in this figure are for 
MDI(A) and SkO' used throughout this work. For con¬ 
trast, we also show results for a typical mean-field the¬ 
oretical model (labelled MFT in the figure) with up to 
quartic scalar self-interactions. The strength parameters 
of this MFT model yield the zero-temperature properties 
of n 0 = 0.155 fm~ 3 , E 0 = 16 MeV, M*/M = 0.7, K 0 = 
222 MeV, S 2 = 30 MeV, and L = 87.0 MeV [HE]. 

In the non-degenerate regime, the exact numerical re¬ 
sults in all cases shown tend to 5/3 as expected. Also 


shown in this figure are results from the expression in 
Eq. (70) which agree very well with the exact results in 


the expected regions of density (that are very nearly in¬ 
dependent of temperature) for all three models (except 
the MFT model at T = 50 MeV, to which we will re¬ 
turn below). Beginning with the non-relativistic models, 
we observe a distinct difference between results for the 
two models in that the MDI(A) model exhibits a pro¬ 
nounced peak for Y p = x = 0.5 whereas the SkO' model 
does not. The origin of these differences can be traced 
back to the behavior of to* ’s of these models with density. 
The presence or absence of a peak can be ascertained by 
examining whether or not 


dr t h dm* 

*r = ° = ^r 


i - 


dm* 
m* dn 


n = 0 (80) 
dn- 


admits a solution. For the MDI(A) model, the solution 
of the above equation occurs at n ~ 0.15 fm -3 at T = 0 
in good agreement with the exact numerical results. As 
the temperature increases toward 50 MeV, the estimate 
from the degenerate limit no longer applies as can be seen 
from the figure. 

We turn now to analyze results of the MFT model, 
particularly in the degenerate region where a peak occurs 
for all values of Y p . In MFT models, the Landau effective 

masses ro* = E* F = \Jp 2 F . + M* 2 , where M* is the Dirac 

effective mass obtained from a self-consistent procedure 
which involves minimizing the energy density (pressure 
at finite T) with respect to density. 

Figure [59] shows M* and ?n* of the MFT model em¬ 
ployed here for T = 0. Note that although M* decreases 
monotonically with density, to* exhibits a minimum and 
rises monotonically with density, a characteristic behav¬ 
ior solely due to relativity (at asymptotic densities, mass 
is overwhelmed by momentum). The density at which 
the minima occur is easily found from the roots of 


Pf dM* 

' -|- 

M* dp Fi 


(81) 


for each Y p . The densities at the minima range from 
0.52-0.57 fm -3 for Y p in the range 0-0.5. This behavior 
of to*’ s is also at the root of the peaks seen in T t h vs n 
as the analysis below shows. 

We will restrict our analysis to the case of PNM as 
the location of the peak in T t .h is not very sensitive to Y p . 
The degenerate limit expressions for the thermal pressure 
and energy density are [551 : 


Pth — 



1 + 



1 - 3 
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M* 



£th = naT 2 , 


(82) 

(83) 


where the level density parameter a = ir 2 E* F /{2p 2 F ). The 
thermal index then becomes 
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FIG. 28. The thermal index Tth vs n for the three models indicated in the figure. Results for the exact calculations are from 
Eqs. ( |49[ ), ( |50| ) and ( |69| ), whereas those for the degenerate limit are from Eqs. ( |70[ )-(|72|. The proton fractions and temperatures 
are indicated in the figures. Results shown are for nucleons only. 



n (fm -3 ) 

FIG. 29. The Dirac and Landau effective masses M* and 
of neutrons in a typical mean-field theoretical (MFT) model 
at zero temperature. 


In the non-relativisitic limit, M*/E* F —► 1 and the log¬ 
arithmic derivative of M* with respect to n tends to 
zero leading to T t h = 5/3. In the ultra-relativistic limit, 
M* —> 0 so that T t h = 4/3. The density at which the 


maximum occurs in T t h can determined from 


dT th _ „ _ _d_ f _p|_ 
dp F dp F V E* 2 


p F M* dM* 
dp F 


E* 2 

a F 


(85) 


The result is n ~ 0.27 fm -3 in good agreement with the 
exact results. Performing the appropriate calculations 
for a mixture of neutrons and protons for Y p in the range 
0.1-0.5 yields densities in the range 0.28-0.31 fm -3 also in 
good agreement with the exact numerical results shown 
Fig. [28} 

The T = 50 MeV results for the MFT model requires 
some explanation. Note that for densities below the peak, 
the exact and degenerate results are somewhat different 
although the qualitative trend is maintained. This is due 
to the fact that in MFT models, M* acquires a tempera¬ 
ture dependence with progressively increasing T, whereas 
the degenerate limit results are calculated with the zero 
temperature M*. 

The preceding analysis sets the stage to assess the role 
leptons (e _ and e + ) play in determining F t ^. As Skyrme 
models lend themselves to a straightforward analysis (re¬ 
sults for other models are similar although quantitative 
differences exist), we show in Fig. 30 'a) the contributions 
from nucleons and leptons to Pth and e t h at T = 5 MeV 
and Y p = 0.5. The ratio Pth/Eth for leptons remains close 
to 1 /3 for all densities with negligible contributions from 
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FIG. 30. (a) Individual contributions from nucleons and lep¬ 
tons to thermal pressure, Pth, and thermal energy density, 
Eth , f° r the SkO' model at the indicated proton fraction and 
temperature, (b) The total Pth and Eth', the inset shows their 
ratio. 


positrons. With increasing density, the ideal gas value of 
Pth/£th = 2/3 for nucleons changes significantly due to 
corrections from the density dependence of m *. From the 
degenerate result for e t h, it is easy to show that its max¬ 
imum value is reached at n = 1 /(2/3) = 0.69 fm -3 using 
ft = 0.724 fm 3 for x = 0.5 for this model. Use of this 
density yields the peak value of e t h — 0.29 MeV fm -3 
in good agreement with the exact result. Likewise, the 
degenerate expression for P t h yields the density at which 
its peak occurs as n = 0.5(1 + \/L8)//3 = 1.62 fm -3 and 
a peak value of P t h — 0.32 MeV fm - , again in good 
agreement with the exact result. 

Figure [30](b) shows the total thermal pressure and en¬ 
ergy density along with its ratio in the inset. It is evident 
that the contributions from the leptons remain subdom¬ 
inant except at very low and very high densities at this 
temperature. With the inclusion of leptons the ratio of 
Pth/£th is less than one at all densities. At large den¬ 
sities, the total thermal pressure approaches a constant 
value whereas £th continues to increase. Thus, the ra¬ 
tio Pt.h/£th, and thereby T t h, decreases at large densities 
resulting in the maximum observed in the inset. The lo¬ 
cation of this maximum can be calculated by solving for 
the density at which 

d f 2,d n Q n T d e Q e 
dn \ 2 a p + a e 

appropriate for SNM with Q n = Q p , a n = a p , and n n = 
n p = n e . Using results from FLT, the above equation 
can be cast in the form 

-? + 7TW^= 0 (87) 

7 dn Q n dn 



with 


= ( 1 Q e m e * \ 
7 V 2 Q n m n *) 



m e * \ 

2 m n * J 


-1 


( 88 ) 


a relation that involves density-dependent effective 
masses of the neutron ( Q n = 5/3 — ?n*/m) and the elec¬ 
tron. A straightforward numerical evaluation yields the 
result ~ 1.37 fm -3 for the density at which the peak in 
T t h occurs in agreement with the exact result. 



FIG. 31. Same as Fig. |30| but at T = 50 MeV. 


For contributions from the leptons ( e~ and e + ), we 
adopt the JEL scheme 96] in all regions of degeneracy 
in our exact numerical calculations. Contributions from 
photons (significant at high values of T and n) are easily 
incorporated using 


7r 2 T 4 
15 (he) 


p _ ^7 

7 3 


and 


4 e 1 

3r’ 


(89) 


respectively. 

In Figs. [TTT]( a.) and (b), results similar to those of 
Fig. [30] are shown, but at T = 50 MeV. Contribu¬ 
tions from the leptons are such that Pth/^t.h remains 
at 1/3 as for T = 5 MeV, with minimal contributions 
from positrons. At very low densities the nucleonic con¬ 
tributions are those of nearly free and non-degenerate 
fermions. With increasing density, however, nucleons en¬ 
ter the semi-degenerate region for which a transparent 
analysis is not possible. While the qualitative features 
are similar to those of Fig. [30] quantitative differences 
are due to the higher temperature in this case. 

Results of T t h for the three models considered in this 
section including the contributions from leptons and pho¬ 
tons are shown in Fig. [32] at the temperatures and pro¬ 
ton fractions shown in the figure. Regions of density and 
temperature for which the degenerate or non-degenerate 
approximation is valid are apparent in this figure. Note 
that the behavior of T t h at supra-nuclear nuclear den¬ 
sities for the non-relativistic models is significantly al¬ 
tered from the results with nucleons only (see Fig. 28 for 
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comparison) primarily because of the contributions from 
leptons. 


B. Adiabatic index 

In hydrodynamics, compressions and rarefactions are 
adiabatic (that is, isentropic) rather than isothermal. 
The adiabatic index defined by 


r s = 


a In P 


5 In n 


n dP 
P dn 


(90) 


is a gauge of the fractional variation of local pressure with 
density and hence the stiffness of the EOS. For our pur¬ 
poses it is more convenient to transform to the variables 
(■ n,T ) which is achieved by the use of Jacobians [51] : 


nfgg d{ppn 

P §fef^(n,T) 

n |f | P d{P,T) 

p m n d M 


r s = 


Cp n dP 

CWP Ihi 


(91) 


The calculation of Tg is facilitated by the isentropes in 
the n-T plane shown in Fig. [33] In contrast to the results 
shown in Fig. |18[ contributions from leptons and photons 
in addition to those from nucleons are included in the 
results of Fig. [33] In order to calculate Tg from Eq. 
(91) one first selects values of {n, T) for which S(n, T) is 
a prescribed constant. These n and T are in turn used 
to calculate P, its derivative (dP/dn)\p, and the specific 
heats Cp and Cy. Alternatively, P can be expressed 
analytically as a function of (S, n) (as for example, in the 
degenerate and non-degenerate limits). We have verified 
that these two approaches yield identical results. 

In the derivation of a degenerate limit expression for 
Tg we must be mindful of the fact that even though 
Cp ~ Cy in FLT, the ratio Cp/Cy cannot be set equal 
to unity. This is because Cp and Cy are polynomials in 
odd powers of T in this limit, and thus 


^ = 1 + C(n)T 2 + ... 


(92) 


where C(n ) is a density dependent factor. The T 2 term 
cannot be ignored at this level of approximation as it 
is of the same order in T as the thermal pressure, P t h- 
We circumvent this problem by turning to Eq. (18), and 
write the adiabatic index as 


r = n ( d(P 0 + P th ) 
S Pq + Pth \ dn 


r 2 CY 


d(P 0 + Pth ) 


dT 


n 


Po + Pth 


9~ 


dPth 


dn 


+ 


T (dP th 


n 2 C v \ dT 


where Po and K are the cold pressure and incompress¬ 
ibility, respectively, the baryonic components of which 
are listed in Appendix |A[ 


Under degenerate conditions (when contributions from 
photons are not significant), Pth can also be expressed in 
terms of S (to leading order) as [551E2] 


Pth — 


nS 2 a iYiQi 

“6“ (Ei a i y i) 


2 , z = n,p, e 


(94) 


which allows Tg to be calculated in a relatively simple 
manner and provides a check on results from Eq. (93). 


When non-degenerate conditions prevail, analytical ex¬ 
pressions for P t h, S, Cy, (dP/dT)\ n , and (dP/dn)\p in 
Eqs. (160), (161), (168), (170) and (171) of Ref. [63] are 
useful, and are employed here. In the semi-degenerate 
region, numerical calculations are unavoidable. 


Results at zero temperature 

As Tg receives significant contributions from the zero 
temperature pressure and its variation with density, we 
begin our analysis by examining Tg = o which highlights 
the role of contributions from electrons. Figure [34] shows 
Tg—o at select values of Y e for the two models employed 
here. The dashed curves in this figure are for nucle¬ 
ons only, whereas the solid curves include the contribu¬ 
tions from leptons. Noteworthy features of the results 
for nucleons are: (i) they diverge at values of n for those 
Y p = Y e for which the pressure vanishes, (ii) they are 
negative in some regions of sub-saturation densities in¬ 
dicating mechanical instability, and (iii) they approach 
constant values for asymptotically low and high densi¬ 
ties. The inclusion of charge-balancing electrons renders 
Tg = o to (i) be positive for all n thereby restoring me¬ 
chanical stability, (ii) exhibit a gradual variation from 
a low to high value with the largest variation occurring 
at sub-saturation densities, and (iii) approach constant 
values for asymptotically low and high densities. These 
features can be quantitatively understood by examining 
the nucleon- and lepton- pressures as functions of density 
(see Fig. 35) as the following analysis shows. 


At zero temperature (entropy), the pressure of rela¬ 
tivistic leptons (electrons) can be written as 


Pi = ^-n 0 (3n 2 n 0 


eye (£) 


4/3 


(95) 


Thus, the adiabatic index of matter with baryons and 
leptons takes the form 


1 2 


r _i ! , n dp b 
S ~° U + Pi dn 


1 + ? 

-n 


-i 


(96) 


Jwhere use of dPi/dn = (4/3)P; has been made. For n not 
too far from no, the pressure of baryons is well approxi- 
(9fk)iated by 


7 Q 

P b ~ Pb(n, x = 1/2) + n 0 (l - 2 x) 2 u 2 —r^- 

du 

P b {n, x = 1/2) ~ u 2 (u — 1), 


(97) 
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FIG. 32. Same as Fig. |28| but with the inclusion of leptons and photons. 



FIG. 33. Isentropes for the MDI(A) and SkO 7 models at the indicated lepton fractions for matter with nucleons, leptons and 
photons. 


where u = n/riQ. As demonstrated below, the above 
relations allow us to gain an analytical understanding of 
the magnitude of T s=o in terms of quantities accessible 
to nuclear experiments. 

For SNM, P b = 0 at u = 1 so that rs =0 takes the 
simple form 


Fs=o(u 


4 (g 0 /9) 

’ 3 + {Pi/no) ' 


(98) 


With Pi = 5.26 MeV frn 3 and values of Kq and no from 


Table. |VI| at u = 1, we obtain rs = o = 2.12(2.09) for the 
MDI(A)[SkO'] model, in very good agreement with the 


exact results in Fig. 34 


The density at which mechanical instability of cold 
baryons-only matter (often referred to as the spinodal 
instability) sets in is determined by dP b /dn = 0. In this 
case, 


f( s p) _ 

1 s=o “ 


4 

3 




-l 


(99) 
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FIG. 34. The zero temperature adiabatic index Ts = o for the MDI(A) and SkO 7 models at at the indicated proton fractions. 
Dashed curves include only nucleons while solid lines include nucleons and leptons. 



FIG. 35. Zero temperature nucleon (dashed curves) and lep¬ 
ton (solid curves) pressures at the indicated proton fractions. 


shows the baryon and lepton pressures vs density from 
which the density regions in which mechanical or spin- 
odal instability occurs for baryons-only matter for the 
two models can be ascertained. For SNM, Eq. (97) im¬ 
plies that M( S p) =2/3, in good agreement with the exact 
results for both models. From Eqs. (95) and (97), the 
leptonic and baryonic pressures at U( sp ) are 


Pi{u (sp )) = 3.06 MeV fm” 3 

Pb(u ( sp )) = -(4/243) K 0 n 0 , (100) 


which yield F^g ~ 1.67(1.65) for the MDI(A)[SkO'] 
model in good accord with the exact results. 


Equation (97) further implies that the density at which 
the derivative of the baryonic pressure, dPb/du is maxi¬ 
mum occurs at u = 1/3. At this density, 


where the superscript ( sp ) stands for spinodal. Figure [35] 


P b (u = 1/3) = -(2/243) I< 0 n 0 
Pi{u = 1/3) = 1.22 MeV fm" 3 . 


( 101 ) 











































32 


The corresponding expression for T s=o is 


r,s=o(w 


1/3) = 


1 K 0 n 0 \ / 2 K 0 n 0 \~ 1 

81 Pi ) 243 Pi ) ’ 


( 102 ) 

which yields a numerical value of ~ 1.28 for both the 
MDI(A) and SkO' models. Note that u = 1/3 marks the 
density at which Ts =0 begins to rise for both models. 
For matter with Y e = 0 (PNM), 


dl\ 

dn 


K(n) 

9 


+ 2 uS ' 2 + u 2 S% , 


(103) 


using which we can write 


Ts=o = 


x 


^(u-l) + S£ 

^ (3u — 2) + 2S ' 2 + uS 2 


(104) 


where the primes above denote derivatives with respect 
to u. The relation above highlights the roles of the first 
and second derivatives of the symmetry energy S 2 (n) at 
sub-saturation densities of SNM. For the special case of 
PNM at u = 1, the adiabatic index therefore becomes 


r S =o(u = l) = 2+\ K °+ Kv . (105) 

O Lj v 

Utilizing the values of Kq, K v and L v in Table [VTl we 
obtain rs =0 (u = 1) = 2.82(2.69) for the MDI(A)[SkO'] 
model in excellent agreement with the exact results of 
2.77 and 2.67, respectively. 

For PNM at u = 2/3, rg =0 simplifies to 


r S =o(« = 2/3)=2(l-^ 


-1 


, 1 S 2 

1 -)-- 

3 S ' 0 


(106) 

where the derivatives above are evaluated at u = 2/3. 
The values of S 2 and S 2 at u = 2/3 are 24.93(-12.15) and 
26.48(-12.9) for the MDI(A) [SkO'], respectively. The 
corresponding values of Ts=o{u = 2/3) are 2.56 and 2.43 
for the two models, which match closely with results of 
exact numerical calculations that yield 2.81 and 2.67, re¬ 
spectively. The differences from the exact results can be 
attributed to terms involving S 4 and its derivatives. 

For values 0 < Y e < 0.5, an analysis similar to that 
presented above can be performed using Eq. (97). The 
density at which the spinodal instability sets in, deter¬ 
mined by solving 


(3„-2) + (1-2x) ! ^ + ^) =0, (107) 

steadily shifts to lower values of u as x decreases toward 
0. Corrections due to to the skewness in E(n,x) and 
contributions from S±(n,x) etc., are small and do not 
affect the qualitative behaviors. As for Y e = 0.5, the 
contribution from electrons to the total pressure entirely 


removes the mechanical (spinodal) instability present in 
baryons-only matter for all Y e as also confirmed by the 
exact numerical results. 

At asymptotically low densities and for Y e ^ 0, rg = o 
approaches 4/3 - the value characteristic for relativistic 
electrons. For PNM, the corresponding value is 5/3, as 
expected for non-relativistic neutrons. For asymptoti¬ 
cally high densities and for all Y e . rs = o is controlled by 
a combination of the highest powers of density in the ex¬ 
pression for the pressure of baryons; both models yield 
the value of ~ 2.5. 


Results for finite entropies 


In Fig. [36j Fg = o and Ts = i for matter with baryons, 
leptons and photons are contrasted for several values of 
Y e . A characteristic feature to note is that Fs = i < 
r s=0 for densities beyond a certain F e -dependent den¬ 
sity, whereas the opposite is the case below that density 
for each Y e . The densities at which rg =1 = Fg =0 occurs 
at sub-saturation densities and ranges from ~ 0.02-0.1 
fm~ 3 for Y e in the range 0-0.5. The approach to the 
corresponding low- and high- density asymptotic values 
also depends on Y e . The largest differences between the 
results of S = 0 and S = 1 are for PNM, the least differ¬ 
ences being those for SNM. 

As S = 1 represents matter that is degenerate for all 
but the very low densities, the ingredients that cause the 
behaviors in Fig. [36] can be identified by utilizing the 
results of FLT. Here we provide relations applicable to 
the cases of SNM and PNM. Those for intermediate Y e ’s 
are straightforward to develop, albeit lengthy. We start 
from 


r s 


n d 
Po + Pth dn 


(Po + Pth) 


(108) 


with Po = Pbo+Pio■ The leptonic pressure and its deriva¬ 
tive at zero temperature (entropy), P;o and dPio/dn , 
are obtained from Eq. (95). The cold baryonic pres¬ 
sure P b0 is that from an appropriate model for baryons 
(MDI(A) or SkO' in this paper), and its derivative 
dPbo/dn = K(n)/9. In the degenerate limit of SNM in 
which p Fb = p Fe = (37r 2 n/2) 1 / 3 , the thermal pressure 
and its derivative are 


Pth 

dPth 

dn 





-2 


Pth 

n 




-1 


d-Qb 2 2 

a b n— - ~a b Q b 

dn 3 



Oe\ _1 " 

4 / ' 

(109) 


The quantities a b = tt 2 ml/(2p^ b ) and a e ~ 7r 2 /(2 p Fc ) 
are the level density parameters of the baryons and rel¬ 
ativistic electrons, respectively. The quantity Q b is as in 
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Eq. (72), and for the SkO' model, 


Q b = ! + l f 1 - — 

2 \ m 


dQb 3 m* b ( ml 

n —— = -—- 1-- 

an 2 mV m 


( 110 ) 

which exemplifies the role of the effective masses in de¬ 
termining P t h■ Also, use of Q e ~ 1/2 has been made in 
obtaining Eq. (109) using which Tg is readily calculated 
for SNM. 

Relations appropriate for PNM are obtained simply by 
setting a e = 0 in Eq. (109) with the result 


Pth — 


S 2 
6 

dPth _ Pth 

dn n 


nQ b 

CLb 


2 _ 
1+ 3 g 


n dQ b 
Q b dn 


( 111 ) 


Consequently, the adiabatic index for PNM becomes 

r s = --- \—+ — •$(-- —)] (112) 

Po + Pth 9 6a 1 3 m ) K 1 


for the SkO' model, which exhibits the role of the ef¬ 
fective mass in determining the thermal component for 
finite S. The result for the MDI(A) model is straight¬ 
forward to obtain, but is somewhat more lengthy due to 
its lengthy expression for the baryon effective mass. We 
have verified that the degenerate limit result for T s ac¬ 
curately reproduces the exact results for all values of Y e 
shown in Fig. |36[ 

Based on a mean-field theoretical model, Shen et al ., 
Ezj have presented results for r,s = i for various Y e 's in¬ 
cluding inhomogeneous phases at sub-saturation densi¬ 
ties, and find that it rises abruptly at the transition from 
non-uniform to uniform matter. In the region of densi¬ 
ties for which an appropriate comparison can be made 
(i.e., the homogeneous phase), our results agree semi- 
quantitatively with those of Ref. m differences being 
due to the different model parameters used. 

Figure [37] shows the effect of increasing S for Y e = 0.5. 
Also shown are results from the degenerate and non¬ 
degenerate limit expressions discussed above. The re¬ 
markable feature of the results shown is that the finite 
and zero entropy adiabatic indices are equal at nearly 
the same density, n x ~ 0.08-0.09 fm~ 3 , for a wide range 
of entropies. For a fixed Y e , the intersection density n x 
can be determined by setting r,s = o = Tg which can be 
rearranged to 


1 dP 0 _ 

P 0 dn \P th 


1 dP t h\ 

dn J 


(113) 


Note that in the degenerate limit, the right hand s ide b e- 
comes independent of S (see the relations in Eq. (109)), 
which explains the near independence with S of n x which 
is mainly determined by the density dependent terms on 
the left hand side and af,,a e and Qb on the right hand 
side. Insofar as the pressure can be expressed as of func¬ 
tion of Yi s and to higher orders in S, a mild dependence 
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FIG. 36. Adiabatic index at fixed entropy at the indicated 
proton fractions for the SkO' and MDI(A) models. Re sults 
at zero and finite entropies are calculated using Eqs. (961 and 
(1081, respectively. 


of n x on S may be expected. Not surprisingly, the non¬ 
degenerate approximation reproduces the exact results 
for Tj vs n better than the degenerate approximation 
with increasing values of S. 

As noted before, the logarithmic derivative 
d In Po/d In n exhibits its largest variation at 
sub-saturation densities (see Fig. 351, whereas 


( d\nP t h/dlnn)\s remains relatively flat (with val¬ 
ues that decrease slightly from ~ 5/3 with increasing 
S ) around the intersection density n x . In addition, 
dlnPo/dlnn > (din P t h/dlnn)\s for n > n x but the 
reverse trend prevails for n < n x . These general features 
for all Y e ’s make Eg < Ts =0 for densities larger than n x 
and Tg > Tg = o for densities lower. 

The degenerate and non-degenerate approximations to 
Tg are compared with the exact numerical results in Fig. 


38 for S = 2, 3, and 4. For both models, and for both 


Ye’s shown, the degenerate approximation reproduces the 
exact results only for S = 2 at supra-nuclear densities. 
As high temperatures are required at these densities for 
S exceeding 2, the non-degenerate approximation fares 
better being indistinguishable from the exact results. 

In Fig. |39[ we show the contributions of nucleons, 
leptons and photons for fixed entropies 5=1 and 4. 
For both proton fractions shown, the dominant contri- 
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FIG. 37. Adiabatic index at fixed entropy from the SkO' and 
MDI(A) models for SNM using Eq. (1081. Results for the 
degenerate (left) and non-degenerate (right) limits are from 
Eqs. (1081 and (1091, respectively. 


butions are from nucleons except at very high densities 
when contributions from leptons begin to become equally 
important. The density-independent contributions from 
photons are significant only for large values of S. 


C. Speed of sound 


We conclude this section with a brief discussion of 
a quantity closely related to Tg, namely the adiabatic 
speed of sound c s given by 



dP 

de 


S,Ye 


(114) 


where c is the speed of light. The energy density e above 
is inclusive of the rest-mass density; in relativistic ap¬ 
proaches this is so by default but in nonrelativistic ones, 
such as in our treatment of nucleons, one must write 
e = e + mn , where the first term here is the internal 
energy density. 


Equation (1141 can be expressed in terms of derivatives 


of the density: 


/ c s N 2 _ dP 


( de 

V c / dn 

S,Ye 

i dn 


s,y b 


dP | 

dn I S,Y e 


E + m + n L v . 

dn IS,Y e 


(115) 


(116) 


Then by taking advantage of the total differential of the 
energy per particle 

P 

dE = TdS -\—=dn + (/ x v + fi e — /i n )dY e (117) 


we transform Eq. (1161 to 


7CsN 2 _ dP 
V c / dn 

= r s - 


S,Ye 

P 


E + m+ — 


- mn 


(118) 

(119) 


thereby making the connection between c s and Tg ex¬ 
plicit. Here, h = e + P is the enthalpy density. Equa¬ 


tion (1191 generalizes the definition of the squared speed 
of sound to finite entropy. This result can also be de¬ 
rived from a time-dependent analysis of traveling sound 
wave s in a me dium (see, e.g, [98, (923). Note that Eqs. 
(1181 and( |119[ ) define c s in terms of previously-calculated 
quantities and, furthermore, allow one to write c s in the 
degenerate limit as an explicit function of the entropy. 

An interesting feature of the degenerate limit is that 
the cold and thermal contributions to c s can be sepa¬ 


rated. We begin from Eq. (1181 which we convert to 


m 


K 

9 


dPth 

dn 


S,Y e 


(H 0 + m) ( 


1 + 


Em t 

Ho+m 


( 120 ) 


where H = E + P/n is the enthalpy per particle. In the 
degenerate region H t h Hq, so that 


m 


K 

9 


dPth I 
dn I S,Y e 


1 - 


Ha+m'j 




Ho + m 


K 

K H t h | 

dPth 1 

9 

9 Ho+m ' 

dn \ S,Y e 




Ho + m 

2 1 

+ yy- 

H 0 + m 


( 121 ) 


( 122 ) 


dPth 

-( C ^) 2 H t d 

dn 

S,Y b ^ C J 


(123) 


In the transition from Eq. (121) to Eq. (122), we have 


discarded the term proportional to dPth/dn x H t h as it 
represents higher-order corrections. 

Results for the squared sound speed vs density for the 
two models are shown in Fig. [40] for values of S up to 
4. Contributions from photons are straightforwardly in¬ 
cluded in Tg, P and h of Eq. (119). Finite entropy contri¬ 


butions to c 2 depend on the density, and are positive up 
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FIG. 38. Adiabatic index at constant entropy for the SkC/ and MDI(A) models using Eqs. (108) and (109). Results for SNM 
[(a) and (c)] and for PNM [(b) and (d)] are shown along with their limiting cases. 
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FIG. 39. Contributions from nucleons, leptons and photons to the total entropy at the indicated proton fractions for the 
MDI(A) and SkO ; models. 


to some high density beyond which the trend is reversed 
as implied in Eq. (122). We note that non-relativistic 
models tend to go acausal at some high density. For a 
prescription to render them causal at finite S , see the 
erratum m to the procedure outlined in Ref. [55] , 


VII. SUMMARY AND CONCLUSION 

In this work, we have sought answers to the following 
two questions, often raised in the community. 

(1) Are equations of state that describe the observed 
collective flow in heavy-ion collisions in the range 0.5- 
2 GeV/A able to account for the largest of the recently 
well-measured neutron star masses ~ 2 Mq? 


(2) How do the thermal properties of such equations of 
state compare with those currently used in simulations 
of core-collapse supernovae, proto-neutron star evolution, 
and binary mergers of compact stars? 

In order to account for the measured mean trans¬ 
verse momentum per nucleon vs rapidity, flow angles, az¬ 
imuthal distributions, and radial flow in medium-energy 
heavy-ion collisions, density and momentum-dependent 
mean fields that matched the real part of the optical 
model potential fits to nucleon-nucleus scattering data 
are required nansnsi. Such mean fields are also the 
result of microscopic many-body calculations of dense 
nucleonic matter [74[ [75]. A characteristic feature of 
these mean fields is that for large momenta, they sat¬ 
urate unlike those given by zero-range Skyrme models 
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FIG. 40. Squared sound speed vs density for various entropy 
per baryon from Eq. (|119[). 


which rise quadratically with momentum (use of Skyrme 
model mean fields yields flow observables that are too 
large compared with data). For rapid use in large-scale 
transport model calculations of heavy-ion collisions, fea¬ 
tures of microscopic calculations have been successfully 
parametrized in schematic models using e.g ., finite-range 
Yukawa-type forces mmm- In this work, we have 
adapted the MDI(0) model of Das et al. [M] in which 
the mean field is both isospin and momentum dependent 
in order to study both the structural aspects of neutron 
stars and the thermal properties of isospin asymmetric 
matter. Our principal findings are summarized below 
along with conclusions. 

The MDI(0) model of Das et al. yields a neutron star 
maximum mass of 1.88 M 0 , which falls slightly below 
the central value of the largest well-measured mass of 
2.01 ± 0.04 M 0 . The radius of the 1.4 M 0 star is 11.27 
km to be compared with 11.5 ± 0.7 km recommended in 
Refs. EM3]. 

To bring the maximum mass closer to the observed 
value, we devised a revised parametrization labeled 
MDI(A) that gave M max = 1.97 M 0 and a radius of 11.58 
km for the 1.4 M 0 star. At the edges of the 1 -a errors 
of the empirical properties A'o, L v and K v at the nu¬ 
clear equilibrium density [with attendant small changes 
in the strength and range parameters of MDI(A)], re¬ 
sults obtained were M. max = 2.15 M 0 and a radius of 
12.13 km for the 1.4 M 0 star. The conclusion from these 
results is that a momentum-dependent EOS that repro¬ 
duces lieavy-ion flow data is well able to yield neutron 
stars in excess of 2 M 0 . 

The pure neutron matter to symmetric nuclear mat¬ 
ter energy difference of the finite range MDI models re¬ 
ceives finite contributions from terms involving potential 
interactions beyond the quadratic term in the neutron- 
proton asymmetry. The magnitudes of such terms are, 
however, relatively small compared to the leading order 
quadratic term. Such beyond-the-leading-order contribu¬ 


tions may also be expected to be present in microscopic 
calculations of isospin asymmetric matter such as BHF, 
DBHF, GFMC, EFT, etc. that employ finite range in¬ 
teractions. The general practice has been to parametrize 
results of these microscopic calculations with a quadratic 
term only. Such potential contributions to Si for l > 2 
are absent in zero-range Skyrme models by construction. 

Our investigations of thermal effects were performed 
for the finite-range MDI(A) model and the zero-range 
Skyrme model called SkO' to provide contrasts. The 
two models give nearly identical results for the energy 
per particle, pressure, chemical potentials, etc., at zero 
temperature for both nuclear matter and pure neutron 
matter at all densities. Consequently, bulk properties 
accessible in experiments involving nuclei near the equi¬ 
librium density of nuclear matter are similar for the two 
models as are their neutron star attributes. However, the 
momentum dependences of their mean field potentials 
are distinctly different, particularly at large momenta. 
The MDI(A) potential saturates logarithmically at high 
momenta whereas the potential of the SkO' model rises 
quadratically with momentum as in all Skyrme models. 
These differences are in turn reflected in their effective 
masses, both in their magnitudes and density depen¬ 
dences. The influence of their differing effective masses 
on the thermal components of energies, pressures, chem¬ 
ical potentials, and the specific heats at constant volume 
and pressure were explored in detail for both symmet¬ 
ric nuclear matter and pure neutron matter. Formulas 
to calculate the thermal components for all proton frac¬ 
tions from 0-0.5 were derived and are presented in the 
text and appendices. 

The isospin splitting of the effective masses in the 
MDI(A) and SkO' models is similar in that to* > to* 
in neutron-rich matter for all densities. Such is not the 
case in many Skyrme models (in which to* < to* for 
neutron-rich matter) which have been successful in pre¬ 
dicting the ground state properties and collective exci¬ 
tations of laboratory nuclei. This reversal in behavior 
of the expected isospin splitting of the neutron and pro¬ 
ton effective masses are exemplified in the finite-range 
MDI(B) model and the zero-range SLy4 model. This 
curious feature led us to establish relations involving the 
strength and range parameters of Skyrme and MDI mod¬ 
els in which to* > to*, a feature exhibited in microscopic 
calculations of dense matter iTHTfSTI . In astrophysical set¬ 
tings where wide ranges of densities, compositions, and 
temperatures are encountered with varying degrees of de¬ 
generacy, the proper disposition of the neutron and pro¬ 
ton effective masses with respect to isospin content is im¬ 
portant insofar they govern all of the thermal variables 
and their response functions. 

The conclusion that emerges from comparing the exact 
numerical results of the finite-range MDI(A) and zero- 
range SkO' models is that the agreement or disagree¬ 
ment depends on the specific thermal variable in ques¬ 
tion as well as the values of all three independent vari¬ 
ables ( n,x,T ). For example, at T = 20 MeV the thermal 
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energies of symmetric nuclear matter begin to differ for 
n > no, whereas they do so for n < 2 n 0 at T=50 MeV. 
For pure neutron matter, there is good agreement for 
all densities at T = 20 MeV, but not so for n > no at 
T = 50 MeV. The thermal pressures on the other hand 
begin to deviate from each other for n > no for both 
x = 0 and a; = 0.5 at T = 20 MeV, whereas they do 
so for n > 2no at T = 50 MeV. Similar dependences 
on ( n,x,T ), albeit in different regimes of (n,x), are seen 
for the thermal chemical potentials and the entropy per 
baryon. The influence of finite range interactions is par¬ 
ticularly evident when one compares the specific heats at 
constant volume, Cy, and at constant pressure, Cp, with 
those of the zero-range interaction model. To varying de¬ 
grees, differences are manifest in Cy for wide ranges of 
n. Likewise, departures in Cp also depend on the values 
of ( n,x,T ). 

Analytical formulas that were derived in the limit¬ 
ing situations of degenerate and non degenerate matter 
shed some light on the origin of the differences described 
above. In the degenerate limit, the differences between 
the magnitudes and density dependences of the nucleon 
effective masses between the two models control the be¬ 
havior of the thermal state variables and their specific 
heats. Due the saturating decline of effective masses with 
density in the MDI models, thermal variables are gener¬ 
ally different in magnitudes than those of the Skyrme 
models. In the non degenerate regime, the influence of 
the finite-range interactions appears to be distinctly dif¬ 
ferent from those of zero-range interactions. 

Some points are worthy of mention. Given the vast 
number of Skyrme parametrizations in the literature, 
many of which are able to successfully predict ground 
state and exitation properties of laboratory nuclei, it is 
quite possible that many of their thermal properties (ex¬ 
ceptions being response functions such as specific heats, 
susceptibilities, etc.) can be made to resemble those of 
models with finite range interactions. However, their in¬ 
ability to describe heavy-ion data satisfactorily would re¬ 
main due to the quadratically rising momentum depen¬ 
dence of their mean field potentials. It must be also be 
noted that at sub-saturation densities, inhomogeneous 
phases of matter must be considered in astrophysical sit¬ 
uations. Whether models with finite-range interactions 
would yield thermal properties of such phases that dif¬ 
fer significantly from those that result from the use of 
zero-range interactions is unexplored and is a subject for 
future study. 

Astrophysical simulations of core-collapse supernovae, 
proto-neutron star evolution and mergers of compact bi¬ 
nary stars are also sensitive to thermal effects in widely 
varying regions of density and temperature. We have 
therefore conducted an in-depth analysis of the thermal 
and adiabatic indices, T t h and Tg (which are measures 
of the stiffness of the equation of state), and the speed of 
sound c s (a guide for the rapidity with which hydrody¬ 
namic evolution proceeds) of hot and dense matter con¬ 
taining nucleons, leptons and photons. Our focus here 


was on the homogeneous phase of matter results of which 
provide a benchmark against which effects of inhomoge¬ 
neous phases known to exist at sub-saturation densities 
and low entropies can be assessed. Studies in this re¬ 
gard are in progress and will be reported separately. Our 
study here has highlighted the influence of nucleon effec¬ 
tive masses, particularly their iso-spin and density depen¬ 
dences, on the thermal properties of baryons. Equally 
important is the role of leptons which serve to remove 
the mechanical (or spinodal) instability of baryons-only 
matter at sub-saturation densities. We find that substan¬ 
tial variations in T t h, Tg and c s begin to occur at sub¬ 
saturation densities before asymptotic values at supra¬ 
nuclear densities are reached. Consequently, laboratory 
experiments and theoretical studies involving neutron- 
rich matter can pin down these quantities. 
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Appendix A: MDI STATE VARIABLES AT T = 0 


In this appendix, analytical expressions for pressure, 
chemical potentials, the nuclear matter incompressibility, 
the symmetry energy and its density derivatives, single 
particle potentials, and the density derivative of the nu¬ 
cleon effective masses at zero temperature are collected. 
We will use the decomposition of the MDI Hamiltonian 
density as H = Hk + Hd, + 'Hm to isolate the contribu¬ 
tions arising from kinetic sources, density-dependent and 
momentum-dependent interactions to the various quan¬ 
tities of interest. 
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1. Pressure P 


3. Nuclear matter incompressibility K 
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4. Symmetry energy S 2 


5. Derivative of the symmetry energy L 2 
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6. Curvature of the symmetry energy K 2 


7. Single-particle potentials Ui 
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8. Derivatives of the effective masses dm*/dm 
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In the last two expressions i ^ j. 


Appendix B: MDI MODELS IN THE 
NON-DEGENERATE LIMIT 


In this appendix we describe the method by which an¬ 
alytical expressions for the state variables valid in the 
limit of small degeneracy (z <C 1) were obtained. 

This method is based primarily on two elements the 
first of which is the assumption that in the nondegenerate 
(ND) limit the system is sufficiently dilute such that the 
term R(p) in the single-particle energy spectrum depends 
weakly on momentum p and therefore the location po of 
the peak of a given thermodynamic integrand is the same 
as that of a free gas: 
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The second element of the method is the happenstance 
that the Taylor series expansion about z = 0 of the state 
functions results in integrals of the form 

I = J dp g(p) (B3) 


which are amenable to estimation via the saddle point 
approximation (10lj : 
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where the primes denote differentiation with respect to p 
and the subscripts 0 evaluation at p = po- The dummy 
variable (3 - to be set equal to 1 at the end of the calcula¬ 
tion - keeps track of the order of the asymptotic expan¬ 
sion. 

From Eq. (B4) it is clear that 
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which is convenient for our purposes because, as will be 
shown subsequently, the integral in Eq. (B6) can be eval¬ 
uated exactly. Henceforth, we shall make use of the def¬ 
inition 
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At this point we must note a caveat with regards to the 
use of the saddle point approximation in the ND limit: 
High temperatures tend to spread the integrands over a 
wider range and thus incur more error in the approxima¬ 
tion. This is not a problem for the temperatures relevant 
to supernovae and compact objects (T < 50 MeV) but 
nevertheless, the method cannot be expected to work for 
arbitrarily high temperatures. Such broadening can often 
be eliminated by a suitable variable change. However, no 
systematic way exists for finding the appropriate trans¬ 
formation. 

In what follows, we demonstrate the method for a 
single-component gas for the sake of simplicity. Gen¬ 
eralization to multi-component gases is straightforward. 
We begin with R(p) given by 
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where C = Ci for pure neutron matter and C = Ci + C u 
for symmetric nuclear matter. Performing the angular 
integrals leads to 
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The assumption of weak dependence on p' allows us to 
treat ln[.. .] 1/,p as constant evaluated a t p 0 R = ( 2 mT ) 1 / 2 
(where only the leading term of Eq. (B2| is kept) and 
take it out of the integral: 
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In the next step, we expand Eq. ( |B11 ) in a Taylor series 
about z = 0 : 
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where terms in the single-particle energy spectrum that 
depend only on the density have been absorbed in z. For 
the first integral we identify 
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where G n i, G n 2 are given by Eq. (B7) and uq = 
(■mT /2 'kH 2 ) 3 G is the quantum concentration. 

By perturbative inversion of Eq. (B18) to second order 
in n/riQ we obtain 
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Thus, here, we effectively keep only the first term in the 
saddle point approximation. This, as well as the trun¬ 
cated Por, induce little error because in the ND limit 
R(p) is itself a small correction to the spectrum. 

Noting that 
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The chemical potential p is related to 2 according to 
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where u(n) are the terms in the single-particle potential 
which depend only on the density. 

We now turn our attention to the kinetic energy den¬ 
sity: 


9 


2 iT 2 h 3 


dp P f P 


(B21) 


1 g 
2 'k 2 H 3 


H(P) 


2H(p) 

dp e T e 


dp e 

— ^ir+4 lnp 


g-jsr+ 41n P 


(B22) 


With the identification 


gn(p) = e R/T , / T i(p) = 7 ^; - 4 hip, p T 1 = (4?nT ) 1/2 

(B23) 

gr 2 (p) = e~ 2R/T , fr 2 {p) = ^-= - 41np, p r2 = (2 mT) 


mT 


we obtain 


r = 3mTG Tl gn Q (z - 


1/2 
(B24) 

(B25) 


= 3 mTn- 


G. 

1 C 

Lx71 


T 1 


1 + 


2 3 / 2 \G n ignQ 


Gn 2 

Gnl 


G 


T 2 


2G t1 , 
(B26) 

The starting point for the calculation of the potential 
energy density is 
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where the factor 2 g would have been g 2 if Ci was equal to 
C u . Invok ing the definition of R(p) [Eq. (B 8 )], we recast 
Eq. (|B27|) as 
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The functions pertaining to the saddle point calculation 
of Gy are 
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Strictly speaking, one should use fy = — In (p 2 f p ) and 


Pv as in Eq. (B2) although at the expense of simplicity. 


With complete expessions for r and V, the total energy 
density is acquired from 
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where Tid is given by Eq. (23). 


The entropy density in this scheme is obtained from 
Eq. (13) which, upon Taylor expansion about z = 0, 
yields 
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Finally, the pressure is given by the thermodynamic iden¬ 
tity 


P = —£ + Ts + pin. 
and the specific heats by 


1 de 
Cv= ndT 


Cp = Cy + 


n 

T 

( OP 

1 dT 

) 2 

n/ 
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dP 


dn 
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(B38) 


We point out that, besides the usual practice of going to 
higher powers of z in n and r, the accuracy of this proce¬ 
dure can be improved by calculating R(p) to 0{P~ l ) in 
the saddle point approximation and by using a beyond- 
leading-order py in the determination of V. 


Appendix C: SPECIFIC HEATS 
Finite-range models 

The main difficulty in calculating Cy and Cp in the 
MDI model arises from the quantity R(p ), which is both 
density- and momentum-dependent. Hence, it cannot be 


absorbed in the chemical potential. To make this state¬ 
ment explicit, consider the spectrum of a single-species 
in the MDI model: 
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2 TO 
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R(p) = 
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The acronyms PNM and SNM above stand for pure neu¬ 
tron mater and symmetric nuclear matter, respectively. 
The last term in Eq. (Cl) has no momentum dependence 


and can thus be subsumed in the chemical potential for 
purposes of integrating over momentum. However, R(p) 
cannot be treated the same way and so the ’’reduced” 
spectrum 
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is implicitly density-dependent causing complications 
when one attempts to take derivatives with respect to 
pi and T because now the Fermi-Dirac distribution has 
the form 
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We obtain dR/dpi\ T and dR/dT | by taking derivatives 

PF 

dp' p' ln 


of Eq. (C2) with respect to the appropriate variables: 
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where Eqs. (C8)-(CIO) have been used for df /dp\ T and 
df/dT | M . 

As dR'/dp\ T and dR'/dT | are momentum- 
dependent, they cannot be t aken o ut of the integrals. We 
must therefore solve Eqs. (Cll )-( [C12 ) self-consistently 
in a manner similar to the one used to compute R itself. 


Equations (|C8|-(C12) are the necessary ingredients for 


calculating the derivatives (with respect to p and T) of 
the kinetic energy density r and the number density n : 
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Here, df/dp\ T and df/dT are given by Eqs. (|C8|) 


(CIOI in which dR/dp\ T and dR/dT\ fi are the self- 
consistent solutions of Eqs. ( |CllD -( [Cl2l ). We also define 
the double integrals (pertaining to the finite-range terms 
of MDI) 
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for the derivatives of which similar considerations hold: 
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Note that a factor of 2 is gained in Eqs. (C19)-(C22) due 
to the interchangeability of / and f. 


For the specific heat at constant volume, we use 
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With the aid of Eqs. (|C4| and (|CT7f , the energy density 

CA 2 r 


of MDI can be expressed as 

e = -—r 4- [ u(n)dn 4- 
2m J 


n 0 


(C25) 


and therefore, its derivatives with respect to p and T are 
given by 
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where Eqs. (C13)-(C16) and (C19)-(C20) are to be uti¬ 


lized. These relations provide a cross-check of evaluat¬ 


ing Cy in Eq. (C24) directly through appropriate tab¬ 


ulations of the non degenerate expression for e (see Ap¬ 
pendix [Bj) as a function of p and T. 

The specific heat at constant pressure is obtained from 


Cp = Cy + 
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where [by means of a Jacobi transformation from (n,T) 
to (p,T)] 
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Starting from Eq. (16), and observing that 
1 dU(n,p) dR 


dp 


dp 
R 

p no 7 t 2 H 3 p 


2C A 2 1 f , , 2 

- / dp p'- 


(A 2 — p 2 +p' 2 ) 


[A 2 + (p + p') 2 ][A 2 4 - (p-p') 2 ] 
we can write the pressure as 


P = -—r 4- [ u(n)dn 4 - 
3m J 


CA 2 T CA 2 r 

h + w— h- 


3 no 


3 n n 


(C31) 


/(C32) 


(C33) 



































































































45 


Hence 
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As with Cy, the derivatives above provide an alter¬ 
nate means to check the evaluation of Cp in Eq. (C28) 


directly through tabulations of the relevant quantities as 
functions of p and T. 


Zero-range models 
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Collecting dn/dp | T and dn/dT\ to the left hand side, 
we get 


For a generic single-species Skyrme model with the 
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With the p and T derivatives of / and n completely de¬ 
termined, the same can be done for r and consequently 
for s, P, and the specific heats. The results are identical 
to those obtained using the method we used in Ref. S3]. 
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